diff options
Diffstat (limited to 'src/math/sqrt.c')
| -rw-r--r-- | src/math/sqrt.c | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/src/math/sqrt.c b/src/math/sqrt.c new file mode 100644 index 00000000..2ebd022b --- /dev/null +++ b/src/math/sqrt.c @@ -0,0 +1,185 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ +/* sqrt(x) + * Return correctly rounded sqrt. + * ------------------------------------------ + * | Use the hardware sqrt if you have one | + * ------------------------------------------ + * Method: + * Bit by bit method using integer arithmetic. (Slow, but portable) + * 1. Normalization + * Scale x to y in [1,4) with even powers of 2: + * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then + * sqrt(x) = 2^k * sqrt(y) + * 2. Bit by bit computation + * Let q = sqrt(y) truncated to i bit after binary point (q = 1), + * i 0 + * i+1 2 + * s = 2*q , and y = 2 * ( y - q ). (1) + * i i i i + * + * To compute q from q , one checks whether + * i+1 i + * + * -(i+1) 2 + * (q + 2 ) <= y. (2) + * i + * -(i+1) + * If (2) is false, then q = q ; otherwise q = q + 2 . + * i+1 i i+1 i + * + * With some algebric manipulation, it is not difficult to see + * that (2) is equivalent to + * -(i+1) + * s + 2 <= y (3) + * i i + * + * The advantage of (3) is that s and y can be computed by + * i i + * the following recurrence formula: + * if (3) is false + * + * s = s , y = y ; (4) + * i+1 i i+1 i + * + * otherwise, + * -i -(i+1) + * s = s + 2 , y = y - s - 2 (5) + * i+1 i i+1 i i + * + * One may easily use induction to prove (4) and (5). + * Note. Since the left hand side of (3) contain only i+2 bits, + * it does not necessary to do a full (53-bit) comparison + * in (3). + * 3. Final rounding + * After generating the 53 bits result, we compute one more bit. + * Together with the remainder, we can decide whether the + * result is exact, bigger than 1/2ulp, or less than 1/2ulp + * (it will never equal to 1/2ulp). + * The rounding mode can be detected by checking whether + * huge + tiny is equal to huge, and whether huge - tiny is + * equal to huge for some floating point number "huge" and "tiny". + * + * Special cases: + * sqrt(+-0) = +-0 ... exact + * sqrt(inf) = inf + * sqrt(-ve) = NaN ... with invalid signal + * sqrt(NaN) = NaN ... with invalid signal for signaling NaN + */ + +#include "libm.h" + +static const double one = 1.0, tiny = 1.0e-300; + +double sqrt(double x) +{ + double z; + int32_t sign = (int)0x80000000; + int32_t ix0,s0,q,m,t,i; + uint32_t r,t1,s1,ix1,q1; + + EXTRACT_WORDS(ix0, ix1, x); + + /* take care of Inf and NaN */ + if ((ix0&0x7ff00000) == 0x7ff00000) { + return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ + } + /* take care of zero */ + if (ix0 <= 0) { + if (((ix0&~sign)|ix1) == 0) + return x; /* sqrt(+-0) = +-0 */ + if (ix0 < 0) + return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ + } + /* normalize x */ + m = ix0>>20; + if (m == 0) { /* subnormal x */ + while (ix0 == 0) { + m -= 21; + ix0 |= (ix1>>11); + ix1 <<= 21; + } + for (i=0; (ix0&0x00100000) == 0; i++) + ix0<<=1; + m -= i - 1; + ix0 |= ix1>>(32-i); + ix1 <<= i; + } + m -= 1023; /* unbias exponent */ + ix0 = (ix0&0x000fffff)|0x00100000; + if (m & 1) { /* odd m, double x to make it even */ + ix0 += ix0 + ((ix1&sign)>>31); + ix1 += ix1; + } + m >>= 1; /* m = [m/2] */ + + /* generate sqrt(x) bit by bit */ + ix0 += ix0 + ((ix1&sign)>>31); + ix1 += ix1; + q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ + r = 0x00200000; /* r = moving bit from right to left */ + + while (r != 0) { + t = s0 + r; + if (t <= ix0) { + s0 = t + r; + ix0 -= t; + q += r; + } + ix0 += ix0 + ((ix1&sign)>>31); + ix1 += ix1; + r >>= 1; + } + + r = sign; + while (r != 0) { + t1 = s1 + r; + t = s0; + if (t < ix0 || (t == ix0 && t1 <= ix1)) { + s1 = t1 + r; + if ((t1&sign) == sign && (s1&sign) == 0) + s0++; + ix0 -= t; + if (ix1 < t1) + ix0--; + ix1 -= t1; + q1 += r; + } + ix0 += ix0 + ((ix1&sign)>>31); + ix1 += ix1; + r >>= 1; + } + + /* use floating add to find out rounding direction */ + if ((ix0|ix1) != 0) { + z = one - tiny; /* raise inexact flag */ + if (z >= one) { + z = one + tiny; + if (q1 == (uint32_t)0xffffffff) { + q1 = 0; + q++; + } else if (z > one) { + if (q1 == (uint32_t)0xfffffffe) + q++; + q1 += 2; + } else + q1 += q1 & 1; + } + } + ix0 = (q>>1) + 0x3fe00000; + ix1 = q1>>1; + if (q&1) + ix1 |= sign; + ix0 += m << 20; + INSERT_WORDS(z, ix0, ix1); + return z; +} |