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 are permitted 6 * provided that this notice is preserved and that due credit is given 7 * to the University of California at Berkeley. The name of the University 8 * may not be used to endorse or promote products derived from this 9 * software without specific prior written permission. This software 10 * is provided ``as is'' without express or implied warranty. 11 * 12 * All recipients should regard themselves as participants in an ongoing 13 * research project and hence should feel obligated to report their 14 * experiences (good or bad) with these elementary function codes, using 15 * the sendbug(8) program, to the authors. 16 */ 17 18 #ifndef lint 19 static char sccsid[] = "@(#)pow.c 5.2 (Berkeley) 04/29/88"; 20 #endif /* not lint */ 21 22 /* POW(X,Y) 23 * RETURN X**Y 24 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 25 * CODED IN C BY K.C. NG, 1/8/85; 26 * REVISED BY K.C. NG on 7/10/85. 27 * 28 * Required system supported functions: 29 * scalb(x,n) 30 * logb(x) 31 * copysign(x,y) 32 * finite(x) 33 * drem(x,y) 34 * 35 * Required kernel functions: 36 * exp__E(a,c) ...return exp(a+c) - 1 - a*a/2 37 * log__L(x) ...return (log(1+x) - 2s)/s, s=x/(2+x) 38 * pow_p(x,y) ...return +(anything)**(finite non zero) 39 * 40 * Method 41 * 1. Compute and return log(x) in three pieces: 42 * log(x) = n*ln2 + hi + lo, 43 * where n is an integer. 44 * 2. Perform y*log(x) by simulating muti-precision arithmetic and 45 * return the answer in three pieces: 46 * y*log(x) = m*ln2 + hi + lo, 47 * where m is an integer. 48 * 3. Return x**y = exp(y*log(x)) 49 * = 2^m * ( exp(hi+lo) ). 50 * 51 * Special cases: 52 * (anything) ** 0 is 1 ; 53 * (anything) ** 1 is itself; 54 * (anything) ** NaN is NaN; 55 * NaN ** (anything except 0) is NaN; 56 * +-(anything > 1) ** +INF is +INF; 57 * +-(anything > 1) ** -INF is +0; 58 * +-(anything < 1) ** +INF is +0; 59 * +-(anything < 1) ** -INF is +INF; 60 * +-1 ** +-INF is NaN and signal INVALID; 61 * +0 ** +(anything except 0, NaN) is +0; 62 * -0 ** +(anything except 0, NaN, odd integer) is +0; 63 * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; 64 * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; 65 * -0 ** (odd integer) = -( +0 ** (odd integer) ); 66 * +INF ** +(anything except 0,NaN) is +INF; 67 * +INF ** -(anything except 0,NaN) is +0; 68 * -INF ** (odd integer) = -( +INF ** (odd integer) ); 69 * -INF ** (even integer) = ( +INF ** (even integer) ); 70 * -INF ** -(anything except integer,NaN) is NaN with signal; 71 * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); 72 * -(anything except 0) ** (non-integer) is NaN with signal; 73 * 74 * Accuracy: 75 * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, 76 * and a Zilog Z8000, 77 * pow(integer,integer) 78 * always returns the correct integer provided it is representable. 79 * In a test run with 100,000 random arguments with 0 < x, y < 20.0 80 * on a VAX, the maximum observed error was 1.79 ulps (units in the 81 * last place). 82 * 83 * Constants : 84 * The hexadecimal values are the intended ones for the following constants. 85 * The decimal values may be used, provided that the compiler will convert 86 * from decimal to binary accurately enough to produce the hexadecimal values 87 * shown. 88 */ 89 90 #if defined(vax)||defined(tahoe) /* VAX D format */ 91 #include <errno.h> 92 extern double infnan(); 93 #ifdef vax 94 #define _0x(A,B) 0x/**/A/**/B 95 #else /* vax */ 96 #define _0x(A,B) 0x/**/B/**/A 97 #endif /* vax */ 98 /* static double */ 99 /* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */ 100 /* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */ 101 /* invln2 = 1.4426950408889634148E0 , Hex 2^ 1 * .B8AA3B295C17F1 */ 102 /* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */ 103 static long ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)}; 104 static long ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)}; 105 static long invln2x[] = { _0x(aa3b,40b8), _0x(17f1,295c)}; 106 static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)}; 107 #define ln2hi (*(double*)ln2hix) 108 #define ln2lo (*(double*)ln2lox) 109 #define invln2 (*(double*)invln2x) 110 #define sqrt2 (*(double*)sqrt2x) 111 #else /* defined(vax)||defined(tahoe) */ 112 static double 113 ln2hi = 6.9314718036912381649E-1 , /*Hex 2^ -1 * 1.62E42FEE00000 */ 114 ln2lo = 1.9082149292705877000E-10 , /*Hex 2^-33 * 1.A39EF35793C76 */ 115 invln2 = 1.4426950408889633870E0 , /*Hex 2^ 0 * 1.71547652B82FE */ 116 sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */ 117 #endif /* defined(vax)||defined(tahoe) */ 118 119 static double zero=0.0, half=1.0/2.0, one=1.0, two=2.0, negone= -1.0; 120 121 double pow(x,y) 122 double x,y; 123 { 124 double drem(),pow_p(),copysign(),t; 125 int finite(); 126 127 if (y==zero) return(one); 128 else if(y==one 129 #if !defined(vax)&&!defined(tahoe) 130 ||x!=x 131 #endif /* !defined(vax)&&!defined(tahoe) */ 132 ) return( x ); /* if x is NaN or y=1 */ 133 #if !defined(vax)&&!defined(tahoe) 134 else if(y!=y) return( y ); /* if y is NaN */ 135 #endif /* !defined(vax)&&!defined(tahoe) */ 136 else if(!finite(y)) /* if y is INF */ 137 if((t=copysign(x,one))==one) return(zero/zero); 138 else if(t>one) return((y>zero)?y:zero); 139 else return((y<zero)?-y:zero); 140 else if(y==two) return(x*x); 141 else if(y==negone) return(one/x); 142 143 /* sign(x) = 1 */ 144 else if(copysign(one,x)==one) return(pow_p(x,y)); 145 146 /* sign(x)= -1 */ 147 /* if y is an even integer */ 148 else if ( (t=drem(y,two)) == zero) return( pow_p(-x,y) ); 149 150 /* if y is an odd integer */ 151 else if (copysign(t,one) == one) return( -pow_p(-x,y) ); 152 153 /* Henceforth y is not an integer */ 154 else if(x==zero) /* x is -0 */ 155 return((y>zero)?-x:one/(-x)); 156 else { /* return NaN */ 157 #if defined(vax)||defined(tahoe) 158 return (infnan(EDOM)); /* NaN */ 159 #else /* defined(vax)||defined(tahoe) */ 160 return(zero/zero); 161 #endif /* defined(vax)||defined(tahoe) */ 162 } 163 } 164 165 /* pow_p(x,y) return x**y for x with sign=1 and finite y */ 166 static double pow_p(x,y) 167 double x,y; 168 { 169 double logb(),scalb(),copysign(),log__L(),exp__E(); 170 double c,s,t,z,tx,ty; 171 #ifdef tahoe 172 double tahoe_tmp; 173 #endif /* tahoe */ 174 float sx,sy; 175 long k=0; 176 int n,m; 177 178 if(x==zero||!finite(x)) { /* if x is +INF or +0 */ 179 #if defined(vax)||defined(tahoe) 180 return((y>zero)?x:infnan(ERANGE)); /* if y<zero, return +INF */ 181 #else /* defined(vax)||defined(tahoe) */ 182 return((y>zero)?x:one/x); 183 #endif /* defined(vax)||defined(tahoe) */ 184 } 185 if(x==1.0) return(x); /* if x=1.0, return 1 since y is finite */ 186 187 /* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */ 188 z=scalb(x,-(n=logb(x))); 189 #if !defined(vax)&&!defined(tahoe) /* IEEE double; subnormal number */ 190 if(n <= -1022) {n += (m=logb(z)); z=scalb(z,-m);} 191 #endif /* !defined(vax)&&!defined(tahoe) */ 192 if(z >= sqrt2 ) {n += 1; z *= half;} z -= one ; 193 194 /* log(x) = nlog2+log(1+z) ~ nlog2 + t + tx */ 195 s=z/(two+z); c=z*z*half; tx=s*(c+log__L(s*s)); 196 t= z-(c-tx); tx += (z-t)-c; 197 198 /* if y*log(x) is neither too big nor too small */ 199 if((s=logb(y)+logb(n+t)) < 12.0) 200 if(s>-60.0) { 201 202 /* compute y*log(x) ~ mlog2 + t + c */ 203 s=y*(n+invln2*t); 204 m=s+copysign(half,s); /* m := nint(y*log(x)) */ 205 k=y; 206 if((double)k==y) { /* if y is an integer */ 207 k = m-k*n; 208 sx=t; tx+=(t-sx); } 209 else { /* if y is not an integer */ 210 k =m; 211 tx+=n*ln2lo; 212 sx=(c=n*ln2hi)+t; tx+=(c-sx)+t; } 213 /* end of checking whether k==y */ 214 215 sy=y; ty=y-sy; /* y ~ sy + ty */ 216 #ifdef tahoe 217 s = (tahoe_tmp = sx)*sy-k*ln2hi; 218 #else /* tahoe */ 219 s=(double)sx*sy-k*ln2hi; /* (sy+ty)*(sx+tx)-kln2 */ 220 #endif /* tahoe */ 221 z=(tx*ty-k*ln2lo); 222 tx=tx*sy; ty=sx*ty; 223 t=ty+z; t+=tx; t+=s; 224 c= -((((t-s)-tx)-ty)-z); 225 226 /* return exp(y*log(x)) */ 227 t += exp__E(t,c); return(scalb(one+t,m)); 228 } 229 /* end of if log(y*log(x)) > -60.0 */ 230 231 else 232 /* exp(+- tiny) = 1 with inexact flag */ 233 {ln2hi+ln2lo; return(one);} 234 else if(copysign(one,y)*(n+invln2*t) <zero) 235 /* exp(-(big#)) underflows to zero */ 236 return(scalb(one,-5000)); 237 else 238 /* exp(+(big#)) overflows to INF */ 239 return(scalb(one, 5000)); 240 241 } 242