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