diff options
Diffstat (limited to 'gcc/real.c')
-rw-r--r-- | gcc/real.c | 78 |
1 files changed, 0 insertions, 78 deletions
@@ -4765,84 +4765,6 @@ const struct real_format real_internal_format = false }; -/* Calculate the square root of X in mode MODE, and store the result - in R. Return TRUE if the operation does not raise an exception. - For details see "High Precision Division and Square Root", - Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June - 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */ - -bool -real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode, - const REAL_VALUE_TYPE *x) -{ - static REAL_VALUE_TYPE halfthree; - static bool init = false; - REAL_VALUE_TYPE h, t, i; - int iter, exp; - - /* sqrt(-0.0) is -0.0. */ - if (real_isnegzero (x)) - { - *r = *x; - return false; - } - - /* Negative arguments return NaN. */ - if (real_isneg (x)) - { - get_canonical_qnan (r, 0); - return false; - } - - /* Infinity and NaN return themselves. */ - if (!real_isfinite (x)) - { - *r = *x; - return false; - } - - if (!init) - { - do_add (&halfthree, &dconst1, &dconsthalf, 0); - init = true; - } - - /* Initial guess for reciprocal sqrt, i. */ - exp = real_exponent (x); - real_ldexp (&i, &dconst1, -exp/2); - - /* Newton's iteration for reciprocal sqrt, i. */ - for (iter = 0; iter < 16; iter++) - { - /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */ - do_multiply (&t, x, &i); - do_multiply (&h, &t, &i); - do_multiply (&t, &h, &dconsthalf); - do_add (&h, &halfthree, &t, 1); - do_multiply (&t, &i, &h); - - /* Check for early convergence. */ - if (iter >= 6 && real_identical (&i, &t)) - break; - - /* ??? Unroll loop to avoid copying. */ - i = t; - } - - /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */ - do_multiply (&t, x, &i); - do_multiply (&h, &t, &i); - do_add (&i, &dconst1, &h, 1); - do_multiply (&h, &t, &i); - do_multiply (&i, &dconsthalf, &h); - do_add (&h, &t, &i, 0); - - /* ??? We need a Tuckerman test to get the last bit. */ - - real_convert (r, mode, &h); - return true; -} - /* Calculate X raised to the integer exponent N in mode MODE and store the result in R. Return true if the result may be inexact due to loss of precision. The algorithm is the classic "left-to-right binary |