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