diff options
| -rw-r--r-- | lib/msun/src/s_cbrt.c | 16 |
1 files changed, 15 insertions, 1 deletions
diff --git a/lib/msun/src/s_cbrt.c b/lib/msun/src/s_cbrt.c index 4e7eaa8..6d90078 100644 --- a/lib/msun/src/s_cbrt.c +++ b/lib/msun/src/s_cbrt.c @@ -72,7 +72,21 @@ cbrt(double x) } else SET_HIGH_WORD(t,sign|(hx/3+B1)); - /* new cbrt to 23 bits; may be implemented in single precision */ + /* + * New cbrt to 26 bits; may be implemented in single precision: + * cbrt(x) = t*cbrt(x/t**3) ~= t*R(x/t**3) + * where R(r) = (14*r**2 + 35*r + 5)/(5*r**2 + 35*r + 14) is the + * (2,2) Pade approximation to cbrt(r) at r = 1. We replace + * r = x/t**3 by 1/r = t**3/x since the latter can be evaluated + * more efficiently, and rearrange the expression for R(r) to use + * 4 additions and 2 divisions instead of the 4 additions, 4 + * multiplications and 1 division that would be required using + * Horner's rule on the numerator and denominator. t being good + * to 32 bits means that |t/cbrt(x)-1| < 1/32, so |x/t**3-1| < 0.1 + * and for R(r) we can use any approximation to cbrt(r) that is good + * to 20 bits on [0.9, 1.1]. The (2,2) Pade approximation is not an + * especially good choice. + */ r=t*t/x; s=C+r*t; t*=G+F/(s+E+D/s); |
