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