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