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