1 /* 2 * Copyright (c) 1985, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 5 * %sccs.include.redist.c% 6 */ 7 8 #ifndef lint 9 static char sccsid[] = "@(#)pow.c 8.1 (Berkeley) 06/04/93"; 10 #endif /* not lint */ 11 12 /* POW(X,Y) 13 * RETURN X**Y 14 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 15 * CODED IN C BY K.C. NG, 1/8/85; 16 * REVISED BY K.C. NG on 7/10/85. 17 * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92. 18 * Required system supported functions: 19 * scalb(x,n) 20 * logb(x) 21 * copysign(x,y) 22 * finite(x) 23 * drem(x,y) 24 * 25 * Required kernel functions: 26 * exp__D(a,c) exp(a + c) for |a| << |c| 27 * struct d_double dlog(x) r.a + r.b, |r.b| < |r.a| 28 * 29 * Method 30 * 1. Compute and return log(x) in three pieces: 31 * log(x) = n*ln2 + hi + lo, 32 * where n is an integer. 33 * 2. Perform y*log(x) by simulating muti-precision arithmetic and 34 * return the answer in three pieces: 35 * y*log(x) = m*ln2 + hi + lo, 36 * where m is an integer. 37 * 3. Return x**y = exp(y*log(x)) 38 * = 2^m * ( exp(hi+lo) ). 39 * 40 * Special cases: 41 * (anything) ** 0 is 1 ; 42 * (anything) ** 1 is itself; 43 * (anything) ** NaN is NaN; 44 * NaN ** (anything except 0) is NaN; 45 * +(anything > 1) ** +INF is +INF; 46 * -(anything > 1) ** +INF is NaN; 47 * +-(anything > 1) ** -INF is +0; 48 * +-(anything < 1) ** +INF is +0; 49 * +(anything < 1) ** -INF is +INF; 50 * -(anything < 1) ** -INF is NaN; 51 * +-1 ** +-INF is NaN and signal INVALID; 52 * +0 ** +(anything except 0, NaN) is +0; 53 * -0 ** +(anything except 0, NaN, odd integer) is +0; 54 * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; 55 * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; 56 * -0 ** (odd integer) = -( +0 ** (odd integer) ); 57 * +INF ** +(anything except 0,NaN) is +INF; 58 * +INF ** -(anything except 0,NaN) is +0; 59 * -INF ** (odd integer) = -( +INF ** (odd integer) ); 60 * -INF ** (even integer) = ( +INF ** (even integer) ); 61 * -INF ** -(anything except integer,NaN) is NaN with signal; 62 * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); 63 * -(anything except 0) ** (non-integer) is NaN with signal; 64 * 65 * Accuracy: 66 * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, 67 * and a Zilog Z8000, 68 * pow(integer,integer) 69 * always returns the correct integer provided it is representable. 70 * In a test run with 100,000 random arguments with 0 < x, y < 20.0 71 * on a VAX, the maximum observed error was 1.79 ulps (units in the 72 * last place). 73 * 74 * Constants : 75 * The hexadecimal values are the intended ones for the following constants. 76 * The decimal values may be used, provided that the compiler will convert 77 * from decimal to binary accurately enough to produce the hexadecimal values 78 * shown. 79 */ 80 81 #include <errno.h> 82 #include <math.h> 83 84 #include "mathimpl.h" 85 86 #if (defined(vax) || defined(tahoe)) 87 #define TRUNC(x) x = (double) (float) x 88 #define _IEEE 0 89 #else 90 #define _IEEE 1 91 #define endian (((*(int *) &one)) ? 1 : 0) 92 #define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000 93 #define infnan(x) 0.0 94 #endif /* vax or tahoe */ 95 96 const static double zero=0.0, one=1.0, two=2.0, negone= -1.0; 97 98 static double pow_P __P((double, double)); 99 100 double pow(x,y) 101 double x,y; 102 { 103 double t; 104 if (y==zero) 105 return (one); 106 else if (y==one || (_IEEE && x != x)) 107 return (x); /* if x is NaN or y=1 */ 108 else if (_IEEE && y!=y) /* if y is NaN */ 109 return (y); 110 else if (!finite(y)) /* if y is INF */ 111 if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */ 112 return (y - y); 113 else if (t>one) 114 return ((y<0)? zero : ((x<zero)? y-y : y)); 115 else 116 return ((y>0)? zero : ((x<0)? y-y : -y)); 117 else if (y==two) 118 return (x*x); 119 else if (y==negone) 120 return (one/x); 121 /* x > 0, x == +0 */ 122 else if (copysign(one, x) == one) 123 return (pow_P(x, y)); 124 125 /* sign(x)= -1 */ 126 /* if y is an even integer */ 127 else if ( (t=drem(y,two)) == zero) 128 return (pow_P(-x, y)); 129 130 /* if y is an odd integer */ 131 else if (copysign(t,one) == one) 132 return (-pow_P(-x, y)); 133 134 /* Henceforth y is not an integer */ 135 else if (x==zero) /* x is -0 */ 136 return ((y>zero)? -x : one/(-x)); 137 else if (_IEEE) 138 return (zero/zero); 139 else 140 return (infnan(EDOM)); 141 } 142 /* kernel function for x >= 0 */ 143 static double 144 #ifdef _ANSI_SOURCE 145 pow_P(double x, double y) 146 #else 147 pow_P(x, y) double x, y; 148 #endif 149 { 150 struct Double s, t, __log__D(); 151 double __exp__D(), huge = 1e300, tiny = 1e-300; 152 153 if (x == zero) 154 if (y > zero) 155 return (zero); 156 else if (_IEEE) 157 return (huge*huge); 158 else 159 return (infnan(ERANGE)); 160 if (x == one) 161 return (one); 162 if (!finite(x)) 163 if (y < zero) 164 return (zero); 165 else if (_IEEE) 166 return (huge*huge); 167 else 168 return (infnan(ERANGE)); 169 if (y >= 7e18) /* infinity */ 170 if (x < 1) 171 return(tiny*tiny); 172 else if (_IEEE) 173 return (huge*huge); 174 else 175 return (infnan(ERANGE)); 176 177 /* Return exp(y*log(x)), using simulated extended */ 178 /* precision for the log and the multiply. */ 179 180 s = __log__D(x); 181 t.a = y; 182 TRUNC(t.a); 183 t.b = y - t.a; 184 t.b = s.b*y + t.b*s.a; 185 t.a *= s.a; 186 s.a = t.a + t.b; 187 s.b = (t.a - s.a) + t.b; 188 return (__exp__D(s.a, s.b)); 189 } 190