/[jscoverage]/trunk/js/src/fdlibm/e_sqrt.c
ViewVC logotype

Annotation of /trunk/js/src/fdlibm/e_sqrt.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2 - (hide annotations)
Wed Aug 1 13:51:53 2007 UTC (14 years, 11 months ago) by siliconforks
File MIME type: text/plain
File size: 16503 byte(s)
Initial import.

1 siliconforks 2 /* -*- 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    

  ViewVC Help
Powered by ViewVC 1.1.24