src/video/e_sqrt.h

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