← Back to branch summary

~ubuntu-branches/ubuntu/precise/kompozer/precise

« back to all changes in this revision

Viewing changes to mozilla/js/src/fdlibm/e_sqrt.c

  • Committer: Bazaar Package Importer
  • Author(s): Anthony Yarusso
  • Date: 2007-08-27 01:11:03 UTC
  • Revision ID: james.westby@ubuntu.com-20070827011103-2jgf4s6532gqu2ka
Tags: upstream-0.7.10
Import upstream version 0.7.10

Show diffs side-by-side

added added

removed removed

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