← Back to branch summary

~ubuntu-branches/ubuntu/precise/classpath/precise

« back to all changes in this revision

Viewing changes to native/fdlibm/e_sqrt.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Koch
  • Date: 2006-05-27 16:11:15 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20060527161115-h6e39eposdt5snb6
Tags: 2:0.91-3
http://bugs.debian.org/368172
* Install header files to /usr/include/classpath.
* debian/control: classpath: Conflict with jamvm < 1.4.3 and
  cacao < 0.96 (Closes: #368172).

Show diffs side-by-side

added added

removed removed

Lines of Context:
native/fdlibm/e_sqrt.c
1
 
 
2
 
/* @(#)e_sqrt.c 5.1 93/09/24 */
 
1
/* @(#)e_sqrt.c 1.3 95/01/18 */
3
2
/*
4
3
 * ====================================================
5
4
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6
5
 *
7
 
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 
6
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8
7
 * Permission to use, copy, modify, and distribute this
9
 
 * software is freely granted, provided that this notice
 
8
 * software is freely granted, provided that this notice 
10
9
 * is preserved.
11
10
 * ====================================================
12
11
 */
16
15
 *           ------------------------------------------
17
16
 *           |  Use the hardware sqrt if you have one |
18
17
 *           ------------------------------------------
19
 
 * Method:
20
 
 *   Bit by bit method using integer arithmetic. (Slow, but portable)
 
18
 * Method: 
 
19
 *   Bit by bit method using integer arithmetic. (Slow, but portable) 
21
20
 *   1. Normalization
22
 
 *      Scale x to y in [1,4) with even powers of 2:
 
21
 *      Scale x to y in [1,4) with even powers of 2: 
23
22
 *      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
24
23
 *              sqrt(x) = 2^k * sqrt(y)
25
24
 *   2. Bit by bit computation
28
27
 *                                     i+1         2
29
28
 *          s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
30
29
 *           i      i            i                 i
31
 
 *
32
 
 *      To compute q    from q , one checks whether
33
 
 *                  i+1       i
 
30
 *                                                        
 
31
 *      To compute q    from q , one checks whether 
 
32
 *                  i+1       i                       
34
33
 *
35
34
 *                            -(i+1) 2
36
35
 *                      (q + 2      ) <= y.                     (2)
40
39
 *                             i+1   i             i+1   i
41
40
 *
42
41
 *      With some algebric manipulation, it is not difficult to see
43
 
 *      that (2) is equivalent to
 
42
 *      that (2) is equivalent to 
44
43
 *                             -(i+1)
45
44
 *                      s  +  2       <= y                      (3)
46
45
 *                       i                i
47
46
 *
48
 
 *      The advantage of (3) is that s  and y  can be computed by
 
47
 *      The advantage of (3) is that s  and y  can be computed by 
49
48
 *                                    i      i
50
49
 *      the following recurrence formula:
51
50
 *          if (3) is false
57
56
 *                         -i                     -(i+1)
58
57
 *          s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
59
58
 *           i+1      i          i+1    i     i
60
 
 *
61
 
 *      One may easily use induction to prove (4) and (5).
 
59
 *                              
 
60
 *      One may easily use induction to prove (4) and (5). 
62
61
 *      Note. Since the left hand side of (3) contain only i+2 bits,
63
 
 *            it does not necessary to do a full (53-bit) comparison
 
62
 *            it does not necessary to do a full (53-bit) comparison 
64
63
 *            in (3).
65
64
 *   3. Final rounding
66
65
 *      After generating the 53 bits result, we compute one more bit.
70
69
 *      The rounding mode can be detected by checking whether
71
70
 *      huge + tiny is equal to huge, and whether huge - tiny is
72
71
 *      equal to huge for some floating point number "huge" and "tiny".
73
 
 *
 
72
 *              
74
73
 * Special cases:
75
74
 *      sqrt(+-0) = +-0         ... exact
76
75
 *      sqrt(inf) = inf
99
98
#endif
100
99
{
101
100
        double z;
102
 
        int32_t sign = (int)0x80000000;
 
101
        int32_t         sign = (int)0x80000000; 
103
102
        uint32_t r,t1,s1,ix1,q1;
104
103
        int32_t ix0,s0,q,m,t,i;
105
104
 
106
105
        EXTRACT_WORDS(ix0,ix1,x);
107
106
 
108
107
    /* take care of Inf and NaN */
109
 
        if((ix0&0x7ff00000)==0x7ff00000) {
 
108
        if((ix0&0x7ff00000)==0x7ff00000) {                      
110
109
            return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
111
110
                                           sqrt(-inf)=sNaN */
112
 
        }
 
111
        } 
113
112
    /* take care of zero */
114
113
        if(ix0<=0) {
115
114
            if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
143
142
        r = 0x00200000;         /* r = moving bit from right to left */
144
143
 
145
144
        while(r!=0) {
146
 
            t = s0+r;
147
 
            if(t<=ix0) {
148
 
                s0   = t+r;
149
 
                ix0 -= t;
150
 
                q   += r;
151
 
            }
 
145
            t = s0+r; 
 
146
            if(t<=ix0) { 
 
147
                s0   = t+r; 
 
148
                ix0 -= t; 
 
149
                q   += r; 
 
150
            } 
152
151
            ix0 += ix0 + ((ix1&sign)>>31);
153
152
            ix1 += ix1;
154
153
            r>>=1;
156
155
 
157
156
        r = sign;
158
157
        while(r!=0) {
159
 
            t1 = s1+r;
 
158
            t1 = s1+r; 
160
159
            t  = s0;
161
 
            if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
 
160
            if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
162
161
                s1  = t1+r;
163
162
                if(((t1&sign)==(uint32_t)sign)&&(s1&sign)==0) s0 += 1;
164
163
                ix0 -= t;
179
178
                if (q1==(uint32_t)0xffffffff) { q1=0; q += 1;}
180
179
                else if (z>one) {
181
180
                    if (q1==(uint32_t)0xfffffffe) q+=1;
182
 
                    q1+=2;
 
181
                    q1+=2; 
183
182
                } else
184
183
                    q1 += (q1&1);
185
184
            }
191
190
        INSERT_WORDS(z,ix0,ix1);
192
191
        return z;
193
192
}
194
 
 
195
193
#endif /* defined(_DOUBLE_IS_32BITS) */
196
194
 
197
195
/*
198
196
Other methods  (use floating-point arithmetic)
199
197
-------------
200
 
(This is a copy of a drafted paper by Prof W. Kahan
 
198
(This is a copy of a drafted paper by Prof W. Kahan 
201
199
and K.C. Ng, written in May, 1986)
202
200
 
203
 
        Two algorithms are given here to implement sqrt(x)
 
201
        Two algorithms are given here to implement sqrt(x) 
204
202
        (IEEE double precision arithmetic) in software.
205
203
        Both supply sqrt(x) correctly rounded. The first algorithm (in
206
204
        Section A) uses newton iterations and involves four divisions.
207
205
        The second one uses reciproot iterations to avoid division, but
208
206
        requires more multiplications. Both algorithms need the ability
209
 
        to chop results of arithmetic operations instead of round them,
 
207
        to chop results of arithmetic operations instead of round them, 
210
208
        and the INEXACT flag to indicate when an arithmetic operation
211
 
        is executed exactly with no roundoff error, all part of the
 
209
        is executed exactly with no roundoff error, all part of the 
212
210
        standard (IEEE 754-1985). The ability to perform shift, add,
213
211
        subtract and logical AND operations upon 32-bit words is needed
214
212
        too, though not part of the standard.
218
216
   (1)  Initial approximation
219
217
 
220
218
        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
 
219
        a floating point number x (in IEEE double format) respectively 
222
220
 
223
221
            1    11                  52                           ...widths
224
222
           ------------------------------------------------------
226
224
           ------------------------------------------------------
227
225
              msb    lsb  msb                                 lsb ...order
228
226
 
229
 
 
 
227
 
230
228
             ------------------------        ------------------------
231
229
        x0:  |s|   e    |    f1     |    x1: |          f2           |
232
230
             ------------------------        ------------------------
251
249
 
252
250
    (2) Iterative refinement
253
251
 
254
 
        Apply Heron's rule three times to y, we have y approximates
 
252
        Apply Heron's rule three times to y, we have y approximates 
255
253
        sqrt(x) to within 1 ulp (Unit in the Last Place):
256
254
 
257
255
                y := (y+x/y)/2          ... almost 17 sig. bits
276
274
        it requires more multiplications and additions. Also x must be
277
275
        scaled in advance to avoid spurious overflow in evaluating the
278
276
        expression 3y*y+x. Hence it is not recommended uless division
279
 
        is slow. If division is very slow, then one should use the
 
277
        is slow. If division is very slow, then one should use the 
280
278
        reciproot algorithm given in section B.
281
279
 
282
280
    (3) Final adjustment
283
281
 
284
 
        By twiddling y's last bit it is possible to force y to be
 
282
        By twiddling y's last bit it is possible to force y to be 
285
283
        correctly rounded according to the prevailing rounding mode
286
284
        as follows. Let r and i be copies of the rounding mode and
287
285
        inexact flag before entering the square root program. Also we
312
310
                I := i;                 ... restore inexact flag
313
311
                R := r;                 ... restore rounded mode
314
312
                return sqrt(x):=y.
315
 
 
 
313
                    
316
314
    (4) Special cases
317
315
 
318
316
        Square root of +inf, +-0, or NaN is itself;
331
329
            k := 0x5fe80000 - (x0>>1);
332
330
            y0:= k - T2[63&(k>>14)].    ... y ~ 1/sqrt(x) to 7.8 bits
333
331
 
334
 
        Here k is a 32-bit integer and T2[] is an integer array
 
332
        Here k is a 32-bit integer and T2[] is an integer array 
335
333
        containing correction terms. Now magically the floating
336
334
        value of y (y's leading 32-bit word is y0, the value of
337
335
        its trailing word y1 is set to zero) approximates 1/sqrt(x)
352
350
 
353
351
        Apply Reciproot iteration three times to y and multiply the
354
352
        result by x to get an approximation z that matches sqrt(x)
355
 
        to about 1 ulp. To be exact, we will have
 
353
        to about 1 ulp. To be exact, we will have 
356
354
                -1ulp < sqrt(x)-z<1.0625ulp.
357
 
 
 
355
        
358
356
        ... set rounding mode to Round-to-nearest
359
357
           y := y*(1.5-0.5*x*y*y)       ... almost 15 sig. bits to 1/sqrt(x)
360
358
           y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
363
361
           z := z + 0.5*z*(1-z*y)       ... about 1 ulp to sqrt(x)
364
362
 
365
363
        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;
 
364
        (a) the term z*y in the final iteration is always less than 1; 
367
365
        (b) the error in the final result is biased upward so that
368
366
                -1 ulp < sqrt(x) - z < 1.0625 ulp
369
367
            instead of |sqrt(x)-z|<1.03125ulp.
370
368
 
371
369
    (3) Final adjustment
372
370
 
373
 
        By twiddling y's last bit it is possible to force y to be
 
371
        By twiddling y's last bit it is possible to force y to be 
374
372
        correctly rounded according to the prevailing rounding mode
375
373
        as follows. Let r and i be copies of the rounding mode and
376
374
        inexact flag before entering the square root program. Also we
410
408
            I := 1;             ... Raise Inexact flag: z is not exact
411
409
        else {
412
410
            j := 1 - [(x0>>20)&1]       ... j = logb(x) mod 2
413
 
            k := z1 >> 26;              ... get z's 25-th and 26-th
 
411
            k := z1 >> 26;              ... get z's 25-th and 26-th 
414
412
                                            fraction bits
415
413
            I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
416
414
        }
417
415
        R:= r           ... restore rounded mode
418
416
        return sqrt(x):=z.
419
417
 
420
 
        If multiplication is cheaper then the foregoing red tape, the
 
418
        If multiplication is cheaper then the foregoing red tape, the 
421
419
        Inexact flag can be evaluated by
422
420
 
423
421
            I := i;
424
422
            I := (z*z!=x) or I.
425
423
 
426
 
        Note that z*z can overwrite I; this value must be sensed if it is
 
424
        Note that z*z can overwrite I; this value must be sensed if it is 
427
425
        True.
428
426
 
429
427
        Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
430
428
        zero.
431
429
 
432
430
                    --------------------
433
 
                z1: |        f2        |
 
431
                z1: |        f2        | 
434
432
                    --------------------
435
433
                bit 31             bit 0
436
434
 
447
445
        11                      01              even
448
446
        -------------------------------------------------
449
447
 
450
 
    (4) Special cases (see (4) of Section A).
451
 
 
 
448
    (4) Special cases (see (4) of Section A).   
 
449
 
452
450
 */
 
451