source:trunk/libs/libmath/e_sqrt.c@469

Last change on this file since 469 was 469, checked in by alain, 3 years ago

1) Introduce the libsemaphore library.
2) Introduce a small libmath library, required by the "fft" application..
3) Introduce the multithreaded "fft" application.
4) Fix a bad synchronisation bug in the Copy-On-Write mechanism.

File size: 13.9 KB
Line
1/*
2 * ====================================================
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/*
13 * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner)
14 */
15
16/* __ieee754_sqrt(x)
17 * Return correctly rounded sqrt.
18 *           ------------------------------------------
19 *           |  Use the hardware sqrt if you have one |
20 *           ------------------------------------------
21 * Method:
22 *   Bit by bit method using integer arithmetic. (Slow, but portable)
23 *   1. Normalization
24 *      Scale x to y in [1,4) with even powers of 2:
25 *      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
26 *              sqrt(x) = 2^k * sqrt(y)
27 *   2. Bit by bit computation
28 *      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
29 *           i                                                   0
30 *                                     i+1         2
31 *          s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
32 *           i      i            i                 i
33 *
34 *      To compute q    from q , one checks whether
35 *                  i+1       i
36 *
37 *                            -(i+1) 2
38 *                      (q + 2      ) <= y.                     (2)
39 *                        i
40 *                                                            -(i+1)
41 *      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
42 *                             i+1   i             i+1   i
43 *
44 *      With some algebric manipulation, it is not difficult to see
45 *      that (2) is equivalent to
46 *                             -(i+1)
47 *                      s  +  2       <= y                      (3)
48 *                       i                i
49 *
50 *      The advantage of (3) is that s  and y  can be computed by
51 *                                    i      i
52 *      the following recurrence formula:
53 *          if (3) is false
54 *
55 *          s     =  s  ,       y    = y   ;                    (4)
56 *           i+1      i          i+1    i
57 *
58 *          otherwise,
59 *                         -i                     -(i+1)
60 *          s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
61 *           i+1      i          i+1    i     i
62 *
63 *      One may easily use induction to prove (4) and (5).
64 *      Note. Since the left hand side of (3) contain only i+2 bits,
65 *            it does not necessary to do a full (53-bit) comparison
66 *            in (3).
67 *   3. Final rounding
68 *      After generating the 53 bits result, we compute one more bit.
69 *      Together with the remainder, we can decide whether the
70 *      result is exact, bigger than 1/2ulp, or less than 1/2ulp
71 *      (it will never equal to 1/2ulp).
72 *      The rounding mode can be detected by checking whether
73 *      huge + tiny is equal to huge, and whether huge - tiny is
74 *      equal to huge for some floating point number "huge" and "tiny".
75 *
76 * Special cases:
77 *      sqrt(+-0) = +-0         ... exact
78 *      sqrt(inf) = inf
79 *      sqrt(-ve) = NaN         ... with invalid signal
80 *      sqrt(NaN) = NaN         ... with invalid signal for signaling NaN
81 *
82 * Other methods : see the appended file at the end of the program below.
83 *---------------
84 */
85
86#include "math.h"
87#include "math_private.h"
88
89static const double one = 1.0, tiny = 1.0e-300;
90
91double __ieee754_sqrt(double x)
92{
93        double z;
94        int32_t sign = (int)0x80000000;
95        int32_t ix0,s0,q,m,t,i;
96        uint32_t r,t1,s1,ix1,q1;
97
98        EXTRACT_WORDS(ix0,ix1,x);
99
100    /* take care of Inf and NaN */
101        if((ix0&0x7ff00000)==0x7ff00000) {
102            return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
103                                           sqrt(-inf)=sNaN */
104        }
105    /* take care of zero */
106        if(ix0<=0) {
107            if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
108            else if(ix0<0)
109                return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
110        }
111    /* normalize x */
112        m = (ix0>>20);
113        if(m==0) {                              /* subnormal x */
114            while(ix0==0) {
115                m -= 21;
116                ix0 |= (ix1>>11); ix1 <<= 21;
117            }
118            for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
119            m -= i-1;
120            ix0 |= (ix1>>(32-i));
121            ix1 <<= i;
122        }
123        m -= 1023;      /* unbias exponent */
124        ix0 = (ix0&0x000fffff)|0x00100000;
125        if(m&1){        /* odd m, double x to make it even */
126            ix0 += ix0 + ((ix1&sign)>>31);
127            ix1 += ix1;
128        }
129        m >>= 1;        /* m = [m/2] */
130
131    /* generate sqrt(x) bit by bit */
132        ix0 += ix0 + ((ix1&sign)>>31);
133        ix1 += ix1;
134        q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
135        r = 0x00200000;         /* r = moving bit from right to left */
136
137        while(r!=0) {
138            t = s0+r;
139            if(t<=ix0) {
140                s0   = t+r;
141                ix0 -= t;
142                q   += r;
143            }
144            ix0 += ix0 + ((ix1&sign)>>31);
145            ix1 += ix1;
146            r>>=1;
147        }
148
149        r = sign;
150        while(r!=0) {
151            t1 = s1+r;
152            t  = s0;
153            if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
154                s1  = t1+r;
155                if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
156                ix0 -= t;
157                if (ix1 < t1) ix0 -= 1;
158                ix1 -= t1;
159                q1  += r;
160            }
161            ix0 += ix0 + ((ix1&sign)>>31);
162            ix1 += ix1;
163            r>>=1;
164        }
165
166    /* use floating add to find out rounding direction */
167        if((ix0|ix1)!=0) {
168            z = one-tiny; /* trigger inexact flag */
169            if (z>=one) {
170                z = one+tiny;
171                if (q1==(uint32_t)0xffffffff) { q1=0; q += 1;}
172                else if (z>one) {
173                    if (q1==(uint32_t)0xfffffffe) q+=1;
174                    q1+=2;
175                } else
176                    q1 += (q1&1);
177            }
178        }
179        ix0 = (q>>1)+0x3fe00000;
180        ix1 =  q1>>1;
181        if ((q&1)==1) ix1 |= sign;
182        ix0 += (m <<20);
183        INSERT_WORDS(z,ix0,ix1);
184        return z;
185}
186
187
188
189/*
190Other methods  (use floating-point arithmetic)
191-------------
192(This is a copy of a drafted paper by Prof W. Kahan
193and 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
208A.  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= {
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
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
314B.  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= {
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,
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
362
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 */
Note: See TracBrowser for help on using the repository browser.