contrib/openbsd_libm/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 /* sqrt(x)
  14  * Return correctly rounded sqrt.
  15  *           ------------------------------------------
  16  *           |  Use the hardware sqrt if you have one |
  17  *           ------------------------------------------
  18  * Method:
  19  *   Bit by bit method using integer arithmetic. (Slow, but portable)
  20  *   1. Normalization
  21  *      Scale x to y in [1,4) with even powers of 2:
  22  *      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
  23  *              sqrt(x) = 2^k * sqrt(y)
  24  *   2. Bit by bit computation
  25  *      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
  26  *           i                                                   0
  27  *                                     i+1         2
  28  *          s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
  29  *           i      i            i                 i
  30  *
  31  *      To compute q    from q , one checks whether
  32  *                  i+1       i
  33  *
  34  *                            -(i+1) 2
  35  *                      (q + 2      ) <= y.                     (2)
  36  *                        i
  37  *                                                            -(i+1)
  38  *      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
  39  *                             i+1   i             i+1   i
  40  *
  41  *      With some algebric manipulation, it is not difficult to see
  42  *      that (2) is equivalent to
  43  *                             -(i+1)
  44  *                      s  +  2       <= y                      (3)
  45  *                       i                i
  46  *
  47  *      The advantage of (3) is that s  and y  can be computed by
  48  *                                    i      i
  49  *      the following recurrence formula:
  50  *          if (3) is false
  51  *
  52  *          s     =  s  ,       y    = y   ;                    (4)
  53  *           i+1      i          i+1    i
  54  *
  55  *          otherwise,
  56  *                         -i                     -(i+1)
  57  *          s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
  58  *           i+1      i          i+1    i     i
  59  *
  60  *      One may easily use induction to prove (4) and (5).
  61  *      Note. Since the left hand side of (3) contain only i+2 bits,
  62  *            it does not necessary to do a full (53-bit) comparison
  63  *            in (3).
  64  *   3. Final rounding
  65  *      After generating the 53 bits result, we compute one more bit.
  66  *      Together with the remainder, we can decide whether the
  67  *      result is exact, bigger than 1/2ulp, or less than 1/2ulp
  68  *      (it will never equal to 1/2ulp).
  69  *      The rounding mode can be detected by checking whether
  70  *      huge + tiny is equal to huge, and whether huge - tiny is
  71  *      equal to huge for some floating point number "huge" and "tiny".
  72  *
  73  * Special cases:
  74  *      sqrt(+-0) = +-0         ... exact
  75  *      sqrt(inf) = inf
  76  *      sqrt(-ve) = NaN         ... with invalid signal
  77  *      sqrt(NaN) = NaN         ... with invalid signal for signaling NaN
  78  *
  79  * Other methods : see the appended file at the end of the program below.
  80  *---------------
  81  */
  82
  83 #include <float.h>
  84 #include <math.h>
  85
  86 #include "math_private.h"
  87
  88 static  const double    one     = 1.0, tiny=1.0e-300;
  89
  90 double
  91 sqrt(double x)
  92 {
  93         double z;
  94         int32_t sign = (int)0x80000000;
  95         int32_t ix0,s0,q,m,t,i;
  96         u_int32_t r,t1,s1,ix1,q1;
  97
  98         EXTRACT_WORDS(ix0,ix1,x);
  99
 100     /* take care of Inf and NaN */
 101         if((ix0&0x7ff00000)==0x7ff00000) {
 102             return x*x+x;               /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
 103                                            sqrt(-inf)=sNaN */
 104         }
 105     /* take care of zero */
 106         if(ix0<=0) {
 107             if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
 108             else if(ix0<0)
 109                 return (x-x)/(x-x);             /* sqrt(-ve) = sNaN */
 110         }
 111     /* normalize x */
 112         m = (ix0>>20);
 113         if(m==0) {                              /* subnormal x */
 114             while(ix0==0) {
 115                 m -= 21;
 116                 ix0 |= (ix1>>11); ix1 <<= 21;
 117             }
 118             for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
 119             m -= i-1;
 120             ix0 |= (ix1>>(32-i));
 121             ix1 <<= i;
 122         }
 123         m -= 1023;      /* unbias exponent */
 124         ix0 = (ix0&0x000fffff)|0x00100000;
 125         if(m&1){        /* odd m, double x to make it even */
 126             ix0 += ix0 + ((ix1&sign)>>31);
 127             ix1 += ix1;
 128         }
 129         m >>= 1;        /* m = [m/2] */
 130
 131     /* generate sqrt(x) bit by bit */
 132         ix0 += ix0 + ((ix1&sign)>>31);
 133         ix1 += ix1;
 134         q = q1 = s0 = s1 = 0;   /* [q,q1] = sqrt(x) */
 135         r = 0x00200000;         /* r = moving bit from right to left */
 136
 137         while(r!=0) {
 138             t = s0+r;
 139             if(t<=ix0) {
 140                 s0   = t+r;
 141                 ix0 -= t;
 142                 q   += r;
 143             }
 144             ix0 += ix0 + ((ix1&sign)>>31);
 145             ix1 += ix1;
 146             r>>=1;
 147         }
 148
 149         r = sign;
 150         while(r!=0) {
 151             t1 = s1+r;
 152             t  = s0;
 153             if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
 154                 s1  = t1+r;
 155                 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
 156                 ix0 -= t;
 157                 if (ix1 < t1) ix0 -= 1;
 158                 ix1 -= t1;
 159                 q1  += r;
 160             }
 161             ix0 += ix0 + ((ix1&sign)>>31);
 162             ix1 += ix1;
 163             r>>=1;
 164         }
 165
 166     /* use floating add to find out rounding direction */
 167         if((ix0|ix1)!=0) {
 168             z = one-tiny; /* trigger inexact flag */
 169             if (z>=one) {
 170                 z = one+tiny;
 171                 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
 172                 else if (z>one) {
 173                     if (q1==(u_int32_t)0xfffffffe) q+=1;
 174                     q1+=2;
 175                 } else
 176                     q1 += (q1&1);
 177             }
 178         }
 179         ix0 = (q>>1)+0x3fe00000;
 180         ix1 =  q1>>1;
 181         if ((q&1)==1) ix1 |= sign;
 182         ix0 += (m <<20);
 183         INSERT_WORDS(z,ix0,ix1);
 184         return z;
 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  */
 443
 444 #if     LDBL_MANT_DIG == DBL_MANT_DIG
 445 __strong_alias(sqrtl, sqrt);
 446 #endif  /* LDBL_MANT_DIG == DBL_MANT_DIG */