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