1 2 /* @(#)e_sqrt.c 1.3 95/01/18 */ 3 /* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 */ 13 14 #include <sys/cdefs.h> 15 __FBSDID("$FreeBSD$"); 16 17 #include <float.h> 18 19 #include "math.h" 20 #include "math_private.h" 21 22 #ifdef USE_BUILTIN_SQRT 23 double 24 __ieee754_sqrt(double x) 25 { 26 return (__builtin_sqrt(x)); 27 } 28 #else 29 /* __ieee754_sqrt(x) 30 * Return correctly rounded sqrt. 31 * ------------------------------------------ 32 * | Use the hardware sqrt if you have one | 33 * ------------------------------------------ 34 * Method: 35 * Bit by bit method using integer arithmetic. (Slow, but portable) 36 * 1. Normalization 37 * Scale x to y in [1,4) with even powers of 2: 38 * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then 39 * sqrt(x) = 2^k * sqrt(y) 40 * 2. Bit by bit computation 41 * Let q = sqrt(y) truncated to i bit after binary point (q = 1), 42 * i 0 43 * i+1 2 44 * s = 2*q , and y = 2 * ( y - q ). (1) 45 * i i i i 46 * 47 * To compute q from q , one checks whether 48 * i+1 i 49 * 50 * -(i+1) 2 51 * (q + 2 ) <= y. (2) 52 * i 53 * -(i+1) 54 * If (2) is false, then q = q ; otherwise q = q + 2 . 55 * i+1 i i+1 i 56 * 57 * With some algebric manipulation, it is not difficult to see 58 * that (2) is equivalent to 59 * -(i+1) 60 * s + 2 <= y (3) 61 * i i 62 * 63 * The advantage of (3) is that s and y can be computed by 64 * i i 65 * the following recurrence formula: 66 * if (3) is false 67 * 68 * s = s , y = y ; (4) 69 * i+1 i i+1 i 70 * 71 * otherwise, 72 * -i -(i+1) 73 * s = s + 2 , y = y - s - 2 (5) 74 * i+1 i i+1 i i 75 * 76 * One may easily use induction to prove (4) and (5). 77 * Note. Since the left hand side of (3) contain only i+2 bits, 78 * it does not necessary to do a full (53-bit) comparison 79 * in (3). 80 * 3. Final rounding 81 * After generating the 53 bits result, we compute one more bit. 82 * Together with the remainder, we can decide whether the 83 * result is exact, bigger than 1/2ulp, or less than 1/2ulp 84 * (it will never equal to 1/2ulp). 85 * The rounding mode can be detected by checking whether 86 * huge + tiny is equal to huge, and whether huge - tiny is 87 * equal to huge for some floating point number "huge" and "tiny". 88 * 89 * Special cases: 90 * sqrt(+-0) = +-0 ... exact 91 * sqrt(inf) = inf 92 * sqrt(-ve) = NaN ... with invalid signal 93 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN 94 * 95 * Other methods : see the appended file at the end of the program below. 96 *--------------- 97 */ 98 99 static const double one = 1.0, tiny=1.0e-300; 100 101 double 102 __ieee754_sqrt(double x) 103 { 104 double z; 105 int32_t sign = (int)0x80000000; 106 int32_t ix0,s0,q,m,t,i; 107 u_int32_t r,t1,s1,ix1,q1; 108 109 EXTRACT_WORDS(ix0,ix1,x); 110 111 /* take care of Inf and NaN */ 112 if((ix0&0x7ff00000)==0x7ff00000) { 113 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 114 sqrt(-inf)=sNaN */ 115 } 116 /* take care of zero */ 117 if(ix0<=0) { 118 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ 119 else if(ix0<0) 120 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ 121 } 122 /* normalize x */ 123 m = (ix0>>20); 124 if(m==0) { /* subnormal x */ 125 while(ix0==0) { 126 m -= 21; 127 ix0 |= (ix1>>11); ix1 <<= 21; 128 } 129 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; 130 m -= i-1; 131 ix0 |= (ix1>>(32-i)); 132 ix1 <<= i; 133 } 134 m -= 1023; /* unbias exponent */ 135 ix0 = (ix0&0x000fffff)|0x00100000; 136 if(m&1){ /* odd m, double x to make it even */ 137 ix0 += ix0 + ((ix1&sign)>>31); 138 ix1 += ix1; 139 } 140 m >>= 1; /* m = [m/2] */ 141 142 /* generate sqrt(x) bit by bit */ 143 ix0 += ix0 + ((ix1&sign)>>31); 144 ix1 += ix1; 145 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ 146 r = 0x00200000; /* r = moving bit from right to left */ 147 148 while(r!=0) { 149 t = s0+r; 150 if(t<=ix0) { 151 s0 = t+r; 152 ix0 -= t; 153 q += r; 154 } 155 ix0 += ix0 + ((ix1&sign)>>31); 156 ix1 += ix1; 157 r>>=1; 158 } 159 160 r = sign; 161 while(r!=0) { 162 t1 = s1+r; 163 t = s0; 164 if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 165 s1 = t1+r; 166 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1; 167 ix0 -= t; 168 if (ix1 < t1) ix0 -= 1; 169 ix1 -= t1; 170 q1 += r; 171 } 172 ix0 += ix0 + ((ix1&sign)>>31); 173 ix1 += ix1; 174 r>>=1; 175 } 176 177 /* use floating add to find out rounding direction */ 178 if((ix0|ix1)!=0) { 179 z = one-tiny; /* trigger inexact flag */ 180 if (z>=one) { 181 z = one+tiny; 182 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;} 183 else if (z>one) { 184 if (q1==(u_int32_t)0xfffffffe) q+=1; 185 q1+=2; 186 } else 187 q1 += (q1&1); 188 } 189 } 190 ix0 = (q>>1)+0x3fe00000; 191 ix1 = q1>>1; 192 if ((q&1)==1) ix1 |= sign; 193 ix0 += (m <<20); 194 INSERT_WORDS(z,ix0,ix1); 195 return z; 196 } 197 #endif 198 199 #if (LDBL_MANT_DIG == 53) 200 __weak_reference(sqrt, sqrtl); 201 #endif 202 203 /* 204 Other methods (use floating-point arithmetic) 205 ------------- 206 (This is a copy of a drafted paper by Prof W. Kahan 207 and K.C. Ng, written in May, 1986) 208 209 Two algorithms are given here to implement sqrt(x) 210 (IEEE double precision arithmetic) in software. 211 Both supply sqrt(x) correctly rounded. The first algorithm (in 212 Section A) uses newton iterations and involves four divisions. 213 The second one uses reciproot iterations to avoid division, but 214 requires more multiplications. Both algorithms need the ability 215 to chop results of arithmetic operations instead of round them, 216 and the INEXACT flag to indicate when an arithmetic operation 217 is executed exactly with no roundoff error, all part of the 218 standard (IEEE 754-1985). The ability to perform shift, add, 219 subtract and logical AND operations upon 32-bit words is needed 220 too, though not part of the standard. 221 222 A. sqrt(x) by Newton Iteration 223 224 (1) Initial approximation 225 226 Let x0 and x1 be the leading and the trailing 32-bit words of 227 a floating point number x (in IEEE double format) respectively 228 229 1 11 52 ...widths 230 ------------------------------------------------------ 231 x: |s| e | f | 232 ------------------------------------------------------ 233 msb lsb msb lsb ...order 234 235 236 ------------------------ ------------------------ 237 x0: |s| e | f1 | x1: | f2 | 238 ------------------------ ------------------------ 239 240 By performing shifts and subtracts on x0 and x1 (both regarded 241 as integers), we obtain an 8-bit approximation of sqrt(x) as 242 follows. 243 244 k := (x0>>1) + 0x1ff80000; 245 y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits 246 Here k is a 32-bit integer and T1[] is an integer array containing 247 correction terms. Now magically the floating value of y (y's 248 leading 32-bit word is y0, the value of its trailing word is 0) 249 approximates sqrt(x) to almost 8-bit. 250 251 Value of T1: 252 static int T1[32]= { 253 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592, 254 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215, 255 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581, 256 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,}; 257 258 (2) Iterative refinement 259 260 Apply Heron's rule three times to y, we have y approximates 261 sqrt(x) to within 1 ulp (Unit in the Last Place): 262 263 y := (y+x/y)/2 ... almost 17 sig. bits 264 y := (y+x/y)/2 ... almost 35 sig. bits 265 y := y-(y-x/y)/2 ... within 1 ulp 266 267 268 Remark 1. 269 Another way to improve y to within 1 ulp is: 270 271 y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x) 272 y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x) 273 274 2 275 (x-y )*y 276 y := y + 2* ---------- ...within 1 ulp 277 2 278 3y + x 279 280 281 This formula has one division fewer than the one above; however, 282 it requires more multiplications and additions. Also x must be 283 scaled in advance to avoid spurious overflow in evaluating the 284 expression 3y*y+x. Hence it is not recommended uless division 285 is slow. If division is very slow, then one should use the 286 reciproot algorithm given in section B. 287 288 (3) Final adjustment 289 290 By twiddling y's last bit it is possible to force y to be 291 correctly rounded according to the prevailing rounding mode 292 as follows. Let r and i be copies of the rounding mode and 293 inexact flag before entering the square root program. Also we 294 use the expression y+-ulp for the next representable floating 295 numbers (up and down) of y. Note that y+-ulp = either fixed 296 point y+-1, or multiply y by nextafter(1,+-inf) in chopped 297 mode. 298 299 I := FALSE; ... reset INEXACT flag I 300 R := RZ; ... set rounding mode to round-toward-zero 301 z := x/y; ... chopped quotient, possibly inexact 302 If(not I) then { ... if the quotient is exact 303 if(z=y) { 304 I := i; ... restore inexact flag 305 R := r; ... restore rounded mode 306 return sqrt(x):=y. 307 } else { 308 z := z - ulp; ... special rounding 309 } 310 } 311 i := TRUE; ... sqrt(x) is inexact 312 If (r=RN) then z=z+ulp ... rounded-to-nearest 313 If (r=RP) then { ... round-toward-+inf 314 y = y+ulp; z=z+ulp; 315 } 316 y := y+z; ... chopped sum 317 y0:=y0-0x00100000; ... y := y/2 is correctly rounded. 318 I := i; ... restore inexact flag 319 R := r; ... restore rounded mode 320 return sqrt(x):=y. 321 322 (4) Special cases 323 324 Square root of +inf, +-0, or NaN is itself; 325 Square root of a negative number is NaN with invalid signal. 326 327 328 B. sqrt(x) by Reciproot Iteration 329 330 (1) Initial approximation 331 332 Let x0 and x1 be the leading and the trailing 32-bit words of 333 a floating point number x (in IEEE double format) respectively 334 (see section A). By performing shifs and subtracts on x0 and y0, 335 we obtain a 7.8-bit approximation of 1/sqrt(x) as follows. 336 337 k := 0x5fe80000 - (x0>>1); 338 y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits 339 340 Here k is a 32-bit integer and T2[] is an integer array 341 containing correction terms. Now magically the floating 342 value of y (y's leading 32-bit word is y0, the value of 343 its trailing word y1 is set to zero) approximates 1/sqrt(x) 344 to almost 7.8-bit. 345 346 Value of T2: 347 static int T2[64]= { 348 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866, 349 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f, 350 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d, 351 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0, 352 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989, 353 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd, 354 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e, 355 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,}; 356 357 (2) Iterative refinement 358 359 Apply Reciproot iteration three times to y and multiply the 360 result by x to get an approximation z that matches sqrt(x) 361 to about 1 ulp. To be exact, we will have 362 -1ulp < sqrt(x)-z<1.0625ulp. 363 364 ... set rounding mode to Round-to-nearest 365 y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x) 366 y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x) 367 ... special arrangement for better accuracy 368 z := x*y ... 29 bits to sqrt(x), with z*y<1 369 z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x) 370 371 Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that 372 (a) the term z*y in the final iteration is always less than 1; 373 (b) the error in the final result is biased upward so that 374 -1 ulp < sqrt(x) - z < 1.0625 ulp 375 instead of |sqrt(x)-z|<1.03125ulp. 376 377 (3) Final adjustment 378 379 By twiddling y's last bit it is possible to force y to be 380 correctly rounded according to the prevailing rounding mode 381 as follows. Let r and i be copies of the rounding mode and 382 inexact flag before entering the square root program. Also we 383 use the expression y+-ulp for the next representable floating 384 numbers (up and down) of y. Note that y+-ulp = either fixed 385 point y+-1, or multiply y by nextafter(1,+-inf) in chopped 386 mode. 387 388 R := RZ; ... set rounding mode to round-toward-zero 389 switch(r) { 390 case RN: ... round-to-nearest 391 if(x<= z*(z-ulp)...chopped) z = z - ulp; else 392 if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp; 393 break; 394 case RZ:case RM: ... round-to-zero or round-to--inf 395 R:=RP; ... reset rounding mod to round-to-+inf 396 if(x<z*z ... rounded up) z = z - ulp; else 397 if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp; 398 break; 399 case RP: ... round-to-+inf 400 if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else 401 if(x>z*z ...chopped) z = z+ulp; 402 break; 403 } 404 405 Remark 3. The above comparisons can be done in fixed point. For 406 example, to compare x and w=z*z chopped, it suffices to compare 407 x1 and w1 (the trailing parts of x and w), regarding them as 408 two's complement integers. 409 410 ...Is z an exact square root? 411 To determine whether z is an exact square root of x, let z1 be the 412 trailing part of z, and also let x0 and x1 be the leading and 413 trailing parts of x. 414 415 If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0 416 I := 1; ... Raise Inexact flag: z is not exact 417 else { 418 j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2 419 k := z1 >> 26; ... get z's 25-th and 26-th 420 fraction bits 421 I := i or (k&j) or ((k&(j+j+1))!=(x1&3)); 422 } 423 R:= r ... restore rounded mode 424 return sqrt(x):=z. 425 426 If multiplication is cheaper then the foregoing red tape, the 427 Inexact flag can be evaluated by 428 429 I := i; 430 I := (z*z!=x) or I. 431 432 Note that z*z can overwrite I; this value must be sensed if it is 433 True. 434 435 Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be 436 zero. 437 438 -------------------- 439 z1: | f2 | 440 -------------------- 441 bit 31 bit 0 442 443 Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd 444 or even of logb(x) have the following relations: 445 446 ------------------------------------------------- 447 bit 27,26 of z1 bit 1,0 of x1 logb(x) 448 ------------------------------------------------- 449 00 00 odd and even 450 01 01 even 451 10 10 odd 452 10 00 even 453 11 01 even 454 ------------------------------------------------- 455 456 (4) Special cases (see (4) of Section A). 457 458 */ 459 460