[388] | 1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
|
---|
| 2 | /*
|
---|
| 3 | * ====================================================
|
---|
| 4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
---|
| 5 | *
|
---|
| 6 | * Developed at SunSoft, a Sun Microsystems, Inc. business.
|
---|
| 7 | * Permission to use, copy, modify, and distribute this
|
---|
| 8 | * software is freely granted, provided that this notice
|
---|
| 9 | * is preserved.
|
---|
| 10 | * ====================================================
|
---|
| 11 | */
|
---|
| 12 | /* sqrt(x)
|
---|
| 13 | * Return correctly rounded sqrt.
|
---|
| 14 | * ------------------------------------------
|
---|
| 15 | * | Use the hardware sqrt if you have one |
|
---|
| 16 | * ------------------------------------------
|
---|
| 17 | * Method:
|
---|
| 18 | * Bit by bit method using integer arithmetic. (Slow, but portable)
|
---|
| 19 | * 1. Normalization
|
---|
| 20 | * Scale x to y in [1,4) with even powers of 2:
|
---|
| 21 | * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
|
---|
| 22 | * sqrt(x) = 2^k * sqrt(y)
|
---|
| 23 | * 2. Bit by bit computation
|
---|
| 24 | * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
|
---|
| 25 | * i 0
|
---|
| 26 | * i+1 2
|
---|
| 27 | * s = 2*q , and y = 2 * ( y - q ). (1)
|
---|
| 28 | * i i i i
|
---|
| 29 | *
|
---|
| 30 | * To compute q from q , one checks whether
|
---|
| 31 | * i+1 i
|
---|
| 32 | *
|
---|
| 33 | * -(i+1) 2
|
---|
| 34 | * (q + 2 ) <= y. (2)
|
---|
| 35 | * i
|
---|
| 36 | * -(i+1)
|
---|
| 37 | * If (2) is false, then q = q ; otherwise q = q + 2 .
|
---|
| 38 | * i+1 i i+1 i
|
---|
| 39 | *
|
---|
| 40 | * With some algebric manipulation, it is not difficult to see
|
---|
| 41 | * that (2) is equivalent to
|
---|
| 42 | * -(i+1)
|
---|
| 43 | * s + 2 <= y (3)
|
---|
| 44 | * i i
|
---|
| 45 | *
|
---|
| 46 | * The advantage of (3) is that s and y can be computed by
|
---|
| 47 | * i i
|
---|
| 48 | * the following recurrence formula:
|
---|
| 49 | * if (3) is false
|
---|
| 50 | *
|
---|
| 51 | * s = s , y = y ; (4)
|
---|
| 52 | * i+1 i i+1 i
|
---|
| 53 | *
|
---|
| 54 | * otherwise,
|
---|
| 55 | * -i -(i+1)
|
---|
| 56 | * s = s + 2 , y = y - s - 2 (5)
|
---|
| 57 | * i+1 i i+1 i i
|
---|
| 58 | *
|
---|
| 59 | * One may easily use induction to prove (4) and (5).
|
---|
| 60 | * Note. Since the left hand side of (3) contain only i+2 bits,
|
---|
| 61 | * it does not necessary to do a full (53-bit) comparison
|
---|
| 62 | * in (3).
|
---|
| 63 | * 3. Final rounding
|
---|
| 64 | * After generating the 53 bits result, we compute one more bit.
|
---|
| 65 | * Together with the remainder, we can decide whether the
|
---|
| 66 | * result is exact, bigger than 1/2ulp, or less than 1/2ulp
|
---|
| 67 | * (it will never equal to 1/2ulp).
|
---|
| 68 | * The rounding mode can be detected by checking whether
|
---|
| 69 | * huge + tiny is equal to huge, and whether huge - tiny is
|
---|
| 70 | * equal to huge for some floating point number "huge" and "tiny".
|
---|
| 71 | *
|
---|
| 72 | * Special cases:
|
---|
| 73 | * sqrt(+-0) = +-0 ... exact
|
---|
| 74 | * sqrt(inf) = inf
|
---|
| 75 | * sqrt(-ve) = NaN ... with invalid signal
|
---|
| 76 | * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
|
---|
| 77 | */
|
---|
| 78 |
|
---|
| 79 | #include "libm.h"
|
---|
| 80 |
|
---|
| 81 | static const double tiny = 1.0e-300;
|
---|
| 82 |
|
---|
| 83 | double sqrt(double x)
|
---|
| 84 | {
|
---|
| 85 | double z;
|
---|
| 86 | int32_t sign = (int)0x80000000;
|
---|
| 87 | int32_t ix0,s0,q,m,t,i;
|
---|
| 88 | uint32_t r,t1,s1,ix1,q1;
|
---|
| 89 |
|
---|
| 90 | EXTRACT_WORDS(ix0, ix1, x);
|
---|
| 91 |
|
---|
| 92 | /* take care of Inf and NaN */
|
---|
| 93 | if ((ix0&0x7ff00000) == 0x7ff00000) {
|
---|
| 94 | return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
|
---|
| 95 | }
|
---|
| 96 | /* take care of zero */
|
---|
| 97 | if (ix0 <= 0) {
|
---|
| 98 | if (((ix0&~sign)|ix1) == 0)
|
---|
| 99 | return x; /* sqrt(+-0) = +-0 */
|
---|
| 100 | if (ix0 < 0)
|
---|
| 101 | return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
|
---|
| 102 | }
|
---|
| 103 | /* normalize x */
|
---|
| 104 | m = ix0>>20;
|
---|
| 105 | if (m == 0) { /* subnormal x */
|
---|
| 106 | while (ix0 == 0) {
|
---|
| 107 | m -= 21;
|
---|
| 108 | ix0 |= (ix1>>11);
|
---|
| 109 | ix1 <<= 21;
|
---|
| 110 | }
|
---|
| 111 | for (i=0; (ix0&0x00100000) == 0; i++)
|
---|
| 112 | ix0<<=1;
|
---|
| 113 | m -= i - 1;
|
---|
| 114 | ix0 |= ix1>>(32-i);
|
---|
| 115 | ix1 <<= i;
|
---|
| 116 | }
|
---|
| 117 | m -= 1023; /* unbias exponent */
|
---|
| 118 | ix0 = (ix0&0x000fffff)|0x00100000;
|
---|
| 119 | if (m & 1) { /* odd m, double x to make it even */
|
---|
| 120 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
| 121 | ix1 += ix1;
|
---|
| 122 | }
|
---|
| 123 | m >>= 1; /* m = [m/2] */
|
---|
| 124 |
|
---|
| 125 | /* generate sqrt(x) bit by bit */
|
---|
| 126 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
| 127 | ix1 += ix1;
|
---|
| 128 | q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
|
---|
| 129 | r = 0x00200000; /* r = moving bit from right to left */
|
---|
| 130 |
|
---|
| 131 | while (r != 0) {
|
---|
| 132 | t = s0 + r;
|
---|
| 133 | if (t <= ix0) {
|
---|
| 134 | s0 = t + r;
|
---|
| 135 | ix0 -= t;
|
---|
| 136 | q += r;
|
---|
| 137 | }
|
---|
| 138 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
| 139 | ix1 += ix1;
|
---|
| 140 | r >>= 1;
|
---|
| 141 | }
|
---|
| 142 |
|
---|
| 143 | r = sign;
|
---|
| 144 | while (r != 0) {
|
---|
| 145 | t1 = s1 + r;
|
---|
| 146 | t = s0;
|
---|
| 147 | if (t < ix0 || (t == ix0 && t1 <= ix1)) {
|
---|
| 148 | s1 = t1 + r;
|
---|
| 149 | if ((t1&sign) == sign && (s1&sign) == 0)
|
---|
| 150 | s0++;
|
---|
| 151 | ix0 -= t;
|
---|
| 152 | if (ix1 < t1)
|
---|
| 153 | ix0--;
|
---|
| 154 | ix1 -= t1;
|
---|
| 155 | q1 += r;
|
---|
| 156 | }
|
---|
| 157 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
| 158 | ix1 += ix1;
|
---|
| 159 | r >>= 1;
|
---|
| 160 | }
|
---|
| 161 |
|
---|
| 162 | /* use floating add to find out rounding direction */
|
---|
| 163 | if ((ix0|ix1) != 0) {
|
---|
| 164 | z = 1.0 - tiny; /* raise inexact flag */
|
---|
| 165 | if (z >= 1.0) {
|
---|
| 166 | z = 1.0 + tiny;
|
---|
| 167 | if (q1 == (uint32_t)0xffffffff) {
|
---|
| 168 | q1 = 0;
|
---|
| 169 | q++;
|
---|
| 170 | } else if (z > 1.0) {
|
---|
| 171 | if (q1 == (uint32_t)0xfffffffe)
|
---|
| 172 | q++;
|
---|
| 173 | q1 += 2;
|
---|
| 174 | } else
|
---|
| 175 | q1 += q1 & 1;
|
---|
| 176 | }
|
---|
| 177 | }
|
---|
| 178 | ix0 = (q>>1) + 0x3fe00000;
|
---|
| 179 | ix1 = q1>>1;
|
---|
| 180 | if (q&1)
|
---|
| 181 | ix1 |= sign;
|
---|
| 182 | ix0 += m << 20;
|
---|
| 183 | INSERT_WORDS(z, ix0, ix1);
|
---|
| 184 | return z;
|
---|
| 185 | }
|
---|