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