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 INSERT_WORDS(z, ix0 + ((uint32_t)m << 20), ix1); 183 return z; 184 } 185