libm/e_sqrt.c

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