StdLib/LibC/Math/e_sqrt.c

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