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