home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / glibc-1.06 / sysdeps / generic / pow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  8.8 KB  |  255 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)pow.c    5.7 (Berkeley) 10/9/90";
  36. #endif /* not lint */
  37.  
  38. /* POW(X,Y)  
  39.  * RETURN X**Y 
  40.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  41.  * CODED IN C BY K.C. NG, 1/8/85; 
  42.  * REVISED BY K.C. NG on 7/10/85.
  43.  *
  44.  * Required system supported functions:
  45.  *      scalb(x,n)      
  46.  *      logb(x)         
  47.  *    copysign(x,y)    
  48.  *    finite(x)    
  49.  *    drem(x,y)
  50.  *
  51.  * Required kernel functions:
  52.  *    exp__E(a,c)    ...return  exp(a+c) - 1 - a*a/2
  53.  *    log__L(x)    ...return  (log(1+x) - 2s)/s, s=x/(2+x) 
  54.  *    pow_p(x,y)    ...return  +(anything)**(finite non zero)
  55.  *
  56.  * Method
  57.  *    1. Compute and return log(x) in three pieces:
  58.  *        log(x) = n*ln2 + hi + lo,
  59.  *       where n is an integer.
  60.  *    2. Perform y*log(x) by simulating muti-precision arithmetic and 
  61.  *       return the answer in three pieces:
  62.  *        y*log(x) = m*ln2 + hi + lo,
  63.  *       where m is an integer.
  64.  *    3. Return x**y = exp(y*log(x))
  65.  *        = 2^m * ( exp(hi+lo) ).
  66.  *
  67.  * Special cases:
  68.  *    (anything) ** 0  is 1 ;
  69.  *    (anything) ** 1  is itself;
  70.  *    (anything) ** NaN is NaN;
  71.  *    NaN ** (anything except 0) is NaN;
  72.  *    +-(anything > 1) ** +INF is +INF;
  73.  *    +-(anything > 1) ** -INF is +0;
  74.  *    +-(anything < 1) ** +INF is +0;
  75.  *    +-(anything < 1) ** -INF is +INF;
  76.  *    +-1 ** +-INF is NaN and signal INVALID;
  77.  *    +0 ** +(anything except 0, NaN)  is +0;
  78.  *    -0 ** +(anything except 0, NaN, odd integer)  is +0;
  79.  *    +0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
  80.  *    -0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
  81.  *    -0 ** (odd integer) = -( +0 ** (odd integer) );
  82.  *    +INF ** +(anything except 0,NaN) is +INF;
  83.  *    +INF ** -(anything except 0,NaN) is +0;
  84.  *    -INF ** (odd integer) = -( +INF ** (odd integer) );
  85.  *    -INF ** (even integer) = ( +INF ** (even integer) );
  86.  *    -INF ** -(anything except integer,NaN) is NaN with signal;
  87.  *    -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
  88.  *    -(anything except 0) ** (non-integer) is NaN with signal;
  89.  *
  90.  * Accuracy:
  91.  *    pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
  92.  *    and a Zilog Z8000,
  93.  *            pow(integer,integer)
  94.  *    always returns the correct integer provided it is representable.
  95.  *    In a test run with 100,000 random arguments with 0 < x, y < 20.0
  96.  *    on a VAX, the maximum observed error was 1.79 ulps (units in the 
  97.  *    last place).
  98.  *
  99.  * Constants :
  100.  * The hexadecimal values are the intended ones for the following constants.
  101.  * The decimal values may be used, provided that the compiler will convert
  102.  * from decimal to binary accurately enough to produce the hexadecimal values
  103.  * shown.
  104.  */
  105.  
  106. #include <errno.h>
  107. #include "mathimpl.h"
  108. #include <limits.h>
  109.  
  110. vc(ln2hi,  6.9314718055829871446E-1  ,7217,4031,0000,f7d0,   0, .B17217F7D00000)
  111. vc(ln2lo,  1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
  112. vc(invln2, 1.4426950408889634148E0   ,aa3b,40b8,17f1,295c,   1, .B8AA3B295C17F1)
  113. vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
  114.  
  115. ic(ln2hi,  6.9314718036912381649E-1,   -1, 1.62E42FEE00000)
  116. ic(ln2lo,  1.9082149292705877000E-10, -33, 1.A39EF35793C76)
  117. ic(invln2, 1.4426950408889633870E0,     0, 1.71547652B82FE)
  118. ic(sqrt2,  1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
  119.  
  120. #ifdef vccast
  121. #define    ln2hi    vccast(ln2hi)
  122. #define    ln2lo    vccast(ln2lo)
  123. #define    invln2    vccast(invln2)
  124. #define    sqrt2    vccast(sqrt2)
  125. #endif
  126.  
  127. static const double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
  128.  
  129. static double pow_p();
  130.  
  131. double pow(x,y)      
  132. double x,y;
  133. {
  134.     double t;
  135.  
  136.     if     (y==zero)      return(one);
  137.     else if(y==one
  138. #if !defined(vax)&&!defined(tahoe)
  139.         || __isnan (x)    /* BSD code did `x!=x' */
  140. #endif    /* !defined(vax)&&!defined(tahoe) */
  141.         ) return( x );      /* if x is NaN or y=1 */
  142. #if !defined(vax)&&!defined(tahoe)
  143.     else if(__isnan(y))         return( y );      /* if y is NaN */
  144. #endif    /* !defined(vax)&&!defined(tahoe) */
  145.     else if(!finite(y))                     /* if y is INF */
  146.          if((t=copysign(x,one))==one) return(zero/zero);
  147.          else if(t>one) return((y>zero)?y:zero);
  148.          else return((y<zero)?-y:zero);
  149.     else if(y==two)       return(x*x);
  150.     else if(y==negone)    return(one/x);
  151.  
  152.     /* sign(x) = 1 */
  153.     else if(copysign(one,x)==one) return(pow_p(x,y));
  154.  
  155.     /* sign(x)= -1 */
  156.     /* if y is an even integer */
  157.     else if ( (t=drem(y,two)) == zero)    return( pow_p(-x,y) );
  158.  
  159.     /* if y is an odd integer */
  160.     else if (copysign(t,one) == one) return( -pow_p(-x,y) );
  161.  
  162.     /* Henceforth y is not an integer */
  163.     else if(x==zero)    /* x is -0 */
  164.         return((y>zero)?-x:one/(-x));
  165.     else {            /* return NaN */
  166. #if defined(vax)||defined(tahoe)
  167.         return (infnan(EDOM));    /* NaN */
  168. #else    /* defined(vax)||defined(tahoe) */
  169.         return(zero/zero);
  170. #endif    /* defined(vax)||defined(tahoe) */
  171.     }
  172. }
  173.  
  174. #ifndef mc68881
  175. /* pow_p(x,y) return x**y for x with sign=1 and finite y */
  176. static double pow_p(x,y)       
  177. double x,y;
  178. {
  179.         double c,s,t,z,tx,ty;
  180. #ifdef tahoe
  181.     double tahoe_tmp;
  182. #endif    /* tahoe */
  183.         float sx,sy;
  184.     long k;
  185.         int n,m;
  186.  
  187.     if(x==zero||!finite(x)) {           /* if x is +INF or +0 */
  188. #if defined(vax)||defined(tahoe)
  189.          return((y>zero)?x:infnan(ERANGE));    /* if y<zero, return +INF */
  190. #else    /* defined(vax)||defined(tahoe) */
  191.          return((y>zero)?x:one/x);
  192. #endif    /* defined(vax)||defined(tahoe) */
  193.     }
  194.     if(x==1.0) return(x);    /* if x=1.0, return 1 since y is finite */
  195.  
  196.     /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
  197.         z=scalb(x,-(n=logb(x)));  
  198. #if !defined(vax)&&!defined(tahoe)    /* IEEE double; subnormal number */
  199.         if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} 
  200. #endif    /* !defined(vax)&&!defined(tahoe) */
  201.         if(z >= sqrt2 ) {n += 1; z *= half;}  z -= one ;
  202.  
  203.     /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
  204.     s=z/(two+z); c=z*z*half; tx=s*(c+log__L(s*s)); 
  205.     t= z-(c-tx); tx += (z-t)-c;
  206.  
  207.    /* if y*log(x) is neither too big nor too small */
  208.     if((s=logb(y)+logb(n+t)) < 12.0) 
  209.         if(s>-60.0) {
  210.  
  211.     /* compute y*log(x) ~ mlog2 + t + c */
  212.             s=y*(n+invln2*t);
  213.                 m=s+copysign(half,s);   /* m := nint(y*log(x)) */ 
  214.         /* (long int) (double) LONG_MIN overflows on some systems.  */
  215.         if (y >= (double) LONG_MIN + 1 && y <= (double) LONG_MAX &&
  216.             (double) (long int) y == y)
  217.           {
  218.             /* Y is an integer */
  219.             k = m - (long int) y * n;
  220.             sx=t; tx+=(t-sx); }
  221.         else    {        /* if y is not an integer */    
  222.             k =m;
  223.              tx+=n*ln2lo;
  224.             sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; }
  225.         /* end of reductions for integral/nonintegral y */
  226.  
  227.                 sy=y; ty=y-sy;          /* y ~ sy + ty */
  228. #ifdef tahoe
  229.         s = (tahoe_tmp = sx)*sy-k*ln2hi;
  230. #else    /* tahoe */
  231.         s=(double)sx*sy-k*ln2hi;        /* (sy+ty)*(sx+tx)-kln2 */
  232. #endif    /* tahoe */
  233.         z=(tx*ty-k*ln2lo);
  234.         tx=tx*sy; ty=sx*ty;
  235.         t=ty+z; t+=tx; t+=s;
  236.         c= -((((t-s)-tx)-ty)-z);
  237.  
  238.         /* return exp(y*log(x)) */
  239.         t += exp__E(t,c); return(scalb(one+t,m));
  240.          }
  241.     /* end of if log(y*log(x)) > -60.0 */
  242.         
  243.         else
  244.         /* exp(+- tiny) = 1 with inexact flag */
  245.             {ln2hi+ln2lo; return(one);}
  246.         else if(copysign(one,y)*(n+invln2*t) <zero)
  247.         /* exp(-(big#)) underflows to zero */
  248.                 return(scalb(one,-5000)); 
  249.         else
  250.             /* exp(+(big#)) overflows to INF */
  251.                 return(scalb(one, 5000)); 
  252.  
  253. }
  254. #endif /* mc68881 */
  255.