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