home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / LIBSRC.ZOO / libsrc / math / pow.c < prev    next >
C/C++ Source or Header  |  1992-02-03  |  9KB  |  251 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.  
  109. vc(ln2hi,  6.9314718055829871446E-1  ,7217,4031,0000,f7d0,   0, .B17217F7D00000)
  110. vc(ln2lo,  1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
  111. vc(invln2, 1.4426950408889634148E0   ,aa3b,40b8,17f1,295c,   1, .B8AA3B295C17F1)
  112. vc(sqrt2,  1.4142135623730950622E0   ,04f3,40b5,de65,33f9,   1, .B504F333F9DE65)
  113.  
  114. ic(ln2hi,  6.9314718036912381649E-1,   -1, 1.62E42FEE00000)
  115. ic(ln2lo,  1.9082149292705877000E-10, -33, 1.A39EF35793C76)
  116. ic(invln2, 1.4426950408889633870E0,     0, 1.71547652B82FE)
  117. ic(sqrt2,  1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
  118.  
  119. #ifdef vccast
  120. #define    ln2hi    vccast(ln2hi)
  121. #define    ln2lo    vccast(ln2lo)
  122. #define    invln2    vccast(invln2)
  123. #define    sqrt2    vccast(sqrt2)
  124. #endif
  125.  
  126. const static double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0;
  127.  
  128. static double pow_p();
  129.  
  130. double pow(x,y)      
  131. double x,y;
  132. {
  133.     double t;
  134.  
  135.     if     (y==zero)      return(one);
  136.     else if(y==one
  137. #if !defined(vax)&&!defined(tahoe)
  138.         ||x!=x
  139. #endif    /* !defined(vax)&&!defined(tahoe) */
  140.         ) return( x );      /* if x is NaN or y=1 */
  141. #if !defined(vax)&&!defined(tahoe)
  142.     else if(y!=y)         return( y );      /* if y is NaN */
  143. #endif    /* !defined(vax)&&!defined(tahoe) */
  144.     else if(!finite(y))                     /* if y is INF */
  145.          if((t=copysign(x,one))==one) return(zero/zero);
  146.          else if(t>one) return((y>zero)?y:zero);
  147.          else return((y<zero)?-y:zero);
  148.     else if(y==two)       return(x*x);
  149.     else if(y==negone)    return(one/x);
  150.  
  151.     /* sign(x) = 1 */
  152.     else if(copysign(one,x)==one) return(pow_p(x,y));
  153.  
  154.     /* sign(x)= -1 */
  155.     /* if y is an even integer */
  156.     else if ( (t=drem(y,two)) == zero)    return( pow_p(-x,y) );
  157.  
  158.     /* if y is an odd integer */
  159.     else if (copysign(t,one) == one) return( -pow_p(-x,y) );
  160.  
  161.     /* Henceforth y is not an integer */
  162.     else if(x==zero)    /* x is -0 */
  163.         return((y>zero)?-x:one/(-x));
  164.     else {            /* return NaN */
  165. #if defined(vax)||defined(tahoe)
  166.         return (infnan(EDOM));    /* NaN */
  167. #else    /* defined(vax)||defined(tahoe) */
  168.         return(zero/zero);
  169. #endif    /* defined(vax)||defined(tahoe) */
  170.     }
  171. }
  172.  
  173. #ifndef mc68881
  174. /* pow_p(x,y) return x**y for x with sign=1 and finite y */
  175. static double pow_p(x,y)       
  176. double x,y;
  177. {
  178.         double c,s,t,z,tx,ty;
  179. #ifdef tahoe
  180.     double tahoe_tmp;
  181. #endif    /* tahoe */
  182.         float sx,sy;
  183.     long k=0;
  184.         int n,m;
  185.  
  186.     if(x==zero||!finite(x)) {           /* if x is +INF or +0 */
  187. #if defined(vax)||defined(tahoe)
  188.          return((y>zero)?x:infnan(ERANGE));    /* if y<zero, return +INF */
  189. #else    /* defined(vax)||defined(tahoe) */
  190.          return((y>zero)?x:one/x);
  191. #endif    /* defined(vax)||defined(tahoe) */
  192.     }
  193.     if(x==1.0) return(x);    /* if x=1.0, return 1 since y is finite */
  194.  
  195.     /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
  196.         z=scalb(x,-(n=logb(x)));  
  197. #if !defined(vax)&&!defined(tahoe)    /* IEEE double; subnormal number */
  198.         if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} 
  199. #endif    /* !defined(vax)&&!defined(tahoe) */
  200.         if(z >= sqrt2 ) {n += 1; z *= half;}  z -= one ;
  201.  
  202.     /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */
  203.     s=z/(two+z); c=z*z*half; tx=s*(c+log__L(s*s)); 
  204.     t= z-(c-tx); tx += (z-t)-c;
  205.  
  206.    /* if y*log(x) is neither too big nor too small */
  207.     if((s=logb(y)+logb(n+t)) < 12.0) 
  208.         if(s>-60.0) {
  209.  
  210.     /* compute y*log(x) ~ mlog2 + t + c */
  211.             s=y*(n+invln2*t);
  212.                 m=s+copysign(half,s);   /* m := nint(y*log(x)) */ 
  213.         k=y; 
  214.         if((double)k==y) {    /* if y is an integer */
  215.             k = m-k*n;
  216.             sx=t; tx+=(t-sx); }
  217.         else    {        /* if y is not an integer */    
  218.             k =m;
  219.              tx+=n*ln2lo;
  220.             sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; }
  221.        /* end of checking whether k==y */
  222.  
  223.                 sy=y; ty=y-sy;          /* y ~ sy + ty */
  224. #ifdef tahoe
  225.         s = (tahoe_tmp = sx)*sy-k*ln2hi;
  226. #else    /* tahoe */
  227.         s=(double)sx*sy-k*ln2hi;        /* (sy+ty)*(sx+tx)-kln2 */
  228. #endif    /* tahoe */
  229.         z=(tx*ty-k*ln2lo);
  230.         tx=tx*sy; ty=sx*ty;
  231.         t=ty+z; t+=tx; t+=s;
  232.         c= -((((t-s)-tx)-ty)-z);
  233.  
  234.         /* return exp(y*log(x)) */
  235.         t += exp__E(t,c); return(scalb(one+t,m));
  236.          }
  237.     /* end of if log(y*log(x)) > -60.0 */
  238.         
  239.         else
  240.         /* exp(+- tiny) = 1 with inexact flag */
  241.             {ln2hi+ln2lo; return(one);}
  242.         else if(copysign(one,y)*(n+invln2*t) <zero)
  243.         /* exp(-(big#)) underflows to zero */
  244.                 return(scalb(one,-5000)); 
  245.         else
  246.             /* exp(+(big#)) overflows to INF */
  247.                 return(scalb(one, 5000)); 
  248.  
  249. }
  250. #endif /* mc68881 */
  251.