gonzui


Format: Advanced Search

t2ex/bsd_source/lib/libc/src_bsd/math/e_sqrt.cbare sourcepermlink (0.01 seconds)

Search this content:

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