/* * Copyright (c) 1985 Regents of the University of California. * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * 3. All advertising materials mentioning features or use of this software * must display the following acknowledgement: * This product includes software developed by the University of * California, Berkeley and its contributors. * 4. Neither the name of the University nor the names of its contributors * may be used to endorse or promote products derived from this software * without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #ifndef lint static char sccsid[] = "@(#)pow.c 5.7 (Berkeley) 10/9/90"; #endif /* not lint */ /* POW(X,Y) * RETURN X**Y * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 7/10/85. * * Required system supported functions: * scalb(x,n) * logb(x) * copysign(x,y) * finite(x) * drem(x,y) * * Required kernel functions: * exp__E(a,c) ...return exp(a+c) - 1 - a*a/2 * log__L(x) ...return (log(1+x) - 2s)/s, s=x/(2+x) * pow_p(x,y) ...return +(anything)**(finite non zero) * * Method * 1. Compute and return log(x) in three pieces: * log(x) = n*ln2 + hi + lo, * where n is an integer. * 2. Perform y*log(x) by simulating muti-precision arithmetic and * return the answer in three pieces: * y*log(x) = m*ln2 + hi + lo, * where m is an integer. * 3. Return x**y = exp(y*log(x)) * = 2^m * ( exp(hi+lo) ). * * Special cases: * (anything) ** 0 is 1 ; * (anything) ** 1 is itself; * (anything) ** NaN is NaN; * NaN ** (anything except 0) is NaN; * +-(anything > 1) ** +INF is +INF; * +-(anything > 1) ** -INF is +0; * +-(anything < 1) ** +INF is +0; * +-(anything < 1) ** -INF is +INF; * +-1 ** +-INF is NaN and signal INVALID; * +0 ** +(anything except 0, NaN) is +0; * -0 ** +(anything except 0, NaN, odd integer) is +0; * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; * -0 ** (odd integer) = -( +0 ** (odd integer) ); * +INF ** +(anything except 0,NaN) is +INF; * +INF ** -(anything except 0,NaN) is +0; * -INF ** (odd integer) = -( +INF ** (odd integer) ); * -INF ** (even integer) = ( +INF ** (even integer) ); * -INF ** -(anything except integer,NaN) is NaN with signal; * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); * -(anything except 0) ** (non-integer) is NaN with signal; * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, * and a Zilog Z8000, * pow(integer,integer) * always returns the correct integer provided it is representable. * In a test run with 100,000 random arguments with 0 < x, y < 20.0 * on a VAX, the maximum observed error was 1.79 ulps (units in the * last place). * * Constants : * The hexadecimal values are the intended ones for the following constants. * The decimal values may be used, provided that the compiler will convert * from decimal to binary accurately enough to produce the hexadecimal values * shown. */ #include #include #include "mathimpl.h" vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000) vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC) vc(invln2, 1.4426950408889634148E0 ,aa3b,40b8,17f1,295c, 1, .B8AA3B295C17F1) vc(sqrt2, 1.4142135623730950622E0 ,04f3,40b5,de65,33f9, 1, .B504F333F9DE65) ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000) ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76) ic(invln2, 1.4426950408889633870E0, 0, 1.71547652B82FE) ic(sqrt2, 1.4142135623730951455E0, 0, 1.6A09E667F3BCD) #ifdef vccast #define ln2hi vccast(ln2hi) #define ln2lo vccast(ln2lo) #define invln2 vccast(invln2) #define sqrt2 vccast(sqrt2) #endif const static double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0; static double pow_p(double, double); double pow(x,y) double x,y; { double t; if (y==zero) return(one); else if(y==one #if !defined(vax)&&!defined(tahoe) ||x!=x #endif /* !defined(vax)&&!defined(tahoe) */ ) return( x ); /* if x is NaN or y=1 */ #if !defined(vax)&&!defined(tahoe) else if(y!=y) return( y ); /* if y is NaN */ #endif /* !defined(vax)&&!defined(tahoe) */ else if(!finite(y)) /* if y is INF */ if((t=copysign(x,one))==one) return(zero/zero); else if(t>one) return((y>zero)?y:zero); else return((yzero)?-x:one/(-x)); else { /* return NaN */ #if defined(vax)||defined(tahoe) return (infnan(EDOM)); /* NaN */ #else /* defined(vax)||defined(tahoe) */ return(zero/zero); #endif /* defined(vax)||defined(tahoe) */ } } #ifndef mc68881 /* pow_p(x,y) return x**y for x with sign=1 and finite y */ static double pow_p(x,y) double x,y; { double c,s,t,z,tx,ty; #ifdef tahoe double tahoe_tmp; #endif /* tahoe */ double errtmp; float sx,sy; long k=0; int n,m; if(x==zero||!finite(x)) { /* if x is +INF or +0 */ #if defined(vax)||defined(tahoe) return((y>zero)?x:infnan(ERANGE)); /* if yzero)?x:one/x); #endif /* defined(vax)||defined(tahoe) */ } if(x==1.0) return(x); /* if x=1.0, return 1 since y is finite */ /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */ z=scalb(x,-(n=logb(x))); #if !defined(vax)&&!defined(tahoe) /* IEEE double; subnormal number */ if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} #endif /* !defined(vax)&&!defined(tahoe) */ if(z >= sqrt2 ) {n += 1; z *= half;} z -= one ; /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */ s=z/(two+z); c=z*z*half; tx=s*(c+log__L(s*s)); t= z-(c-tx); tx += (z-t)-c; /* if y*log(x) is neither too big nor too small */ if((s=logb(y)+logb(n+t)) < 12.0) if(s>-60.0) { /* compute y*log(x) ~ mlog2 + t + c */ s=y*(n+invln2*t); m=s+copysign(half,s); /* m := nint(y*log(x)) */ k=y; if(y > (double)LONG_MIN && y < (double)LONG_MAX && (double)(long)y==y) { /* y is an integer */ k = m-(long)y*n; sx=t; tx+=(t-sx); } else { /* if y is not an integer */ k =m; tx+=n*ln2lo; sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; } /* end of checking whether k==y */ sy=y; ty=y-sy; /* y ~ sy + ty */ #ifdef tahoe s = (tahoe_tmp = sx)*sy-k*ln2hi; #else /* tahoe */ s=(double)sx*sy-k*ln2hi; /* (sy+ty)*(sx+tx)-kln2 */ #endif /* tahoe */ z=(tx*ty-k*ln2lo); tx=tx*sy; ty=sx*ty; t=ty+z; t+=tx; t+=s; c= -((((t-s)-tx)-ty)-z); /* return exp(y*log(x)) */ t += exp__E(t,c); return(scalb(one+t,m)); } /* end of if log(y*log(x)) > -60.0 */ else /* exp(+- tiny) = 1 with inexact flag */ {errtmp=ln2hi+ln2lo; return(one);} else if(copysign(one,y)*(n+invln2*t)