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