e_sqrt.c - vx32 - Local 9vx git repository for patches.
(HTM) git clone git://r-36.net/vx32
(DIR) Log
(DIR) Files
(DIR) Refs
---
e_sqrt.c (14308B)
---
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 #ifndef lint
14 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_sqrt.c,v 1.9 2003/07/23 04:53:46 peter 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 static const double one = 1.0, tiny=1.0e-300;
91
92 double
93 __ieee754_sqrt(double x)
94 {
95 double z;
96 int32_t sign = (int)0x80000000;
97 int32_t ix0,s0,q,m,t,i;
98 u_int32_t r,t1,s1,ix1,q1;
99
100 EXTRACT_WORDS(ix0,ix1,x);
101
102 /* take care of Inf and NaN */
103 if((ix0&0x7ff00000)==0x7ff00000) {
104 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
105 sqrt(-inf)=sNaN */
106 }
107 /* take care of zero */
108 if(ix0<=0) {
109 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
110 else if(ix0<0)
111 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
112 }
113 /* normalize x */
114 m = (ix0>>20);
115 if(m==0) { /* subnormal x */
116 while(ix0==0) {
117 m -= 21;
118 ix0 |= (ix1>>11); ix1 <<= 21;
119 }
120 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
121 m -= i-1;
122 ix0 |= (ix1>>(32-i));
123 ix1 <<= i;
124 }
125 m -= 1023; /* unbias exponent */
126 ix0 = (ix0&0x000fffff)|0x00100000;
127 if(m&1){ /* odd m, double x to make it even */
128 ix0 += ix0 + ((ix1&sign)>>31);
129 ix1 += ix1;
130 }
131 m >>= 1; /* m = [m/2] */
132
133 /* generate sqrt(x) bit by bit */
134 ix0 += ix0 + ((ix1&sign)>>31);
135 ix1 += ix1;
136 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
137 r = 0x00200000; /* r = moving bit from right to left */
138
139 while(r!=0) {
140 t = s0+r;
141 if(t<=ix0) {
142 s0 = t+r;
143 ix0 -= t;
144 q += r;
145 }
146 ix0 += ix0 + ((ix1&sign)>>31);
147 ix1 += ix1;
148 r>>=1;
149 }
150
151 r = sign;
152 while(r!=0) {
153 t1 = s1+r;
154 t = s0;
155 if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
156 s1 = t1+r;
157 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
158 ix0 -= t;
159 if (ix1 < t1) ix0 -= 1;
160 ix1 -= t1;
161 q1 += r;
162 }
163 ix0 += ix0 + ((ix1&sign)>>31);
164 ix1 += ix1;
165 r>>=1;
166 }
167
168 /* use floating add to find out rounding direction */
169 if((ix0|ix1)!=0) {
170 z = one-tiny; /* trigger inexact flag */
171 if (z>=one) {
172 z = one+tiny;
173 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
174 else if (z>one) {
175 if (q1==(u_int32_t)0xfffffffe) q+=1;
176 q1+=2;
177 } else
178 q1 += (q1&1);
179 }
180 }
181 ix0 = (q>>1)+0x3fe00000;
182 ix1 = q1>>1;
183 if ((q&1)==1) ix1 |= sign;
184 ix0 += (m <<20);
185 INSERT_WORDS(z,ix0,ix1);
186 return z;
187 }
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 uless 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