lib/msun/src/e_sqrt.c

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