1 /* $NetBSD: n_pow.c,v 1.6 2002/06/15 00:10:18 matt Exp $ */ 2 /* 3 * Copyright (c) 1985, 1993 4 * The Regents of the University of California. All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 3. All advertising materials mentioning features or use of this software 15 * must display the following acknowledgement: 16 * This product includes software developed by the University of 17 * California, Berkeley and its contributors. 18 * 4. Neither the name of the University nor the names of its contributors 19 * may be used to endorse or promote products derived from this software 20 * without specific prior written permission. 21 * 22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32 * SUCH DAMAGE. 33 */ 34 35 #ifndef lint 36 #if 0 37 static char sccsid[] = "@(#)pow.c 8.1 (Berkeley) 6/4/93"; 38 #endif 39 #endif /* not lint */ 40 41 /* POW(X,Y) 42 * RETURN X**Y 43 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 44 * CODED IN C BY K.C. NG, 1/8/85; 45 * REVISED BY K.C. NG on 7/10/85. 46 * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92. 47 * Required system supported functions: 48 * scalb(x,n) 49 * logb(x) 50 * copysign(x,y) 51 * finite(x) 52 * drem(x,y) 53 * 54 * Required kernel functions: 55 * exp__D(a,c) exp(a + c) for |a| << |c| 56 * struct d_double dlog(x) r.a + r.b, |r.b| < |r.a| 57 * 58 * Method 59 * 1. Compute and return log(x) in three pieces: 60 * log(x) = n*ln2 + hi + lo, 61 * where n is an integer. 62 * 2. Perform y*log(x) by simulating muti-precision arithmetic and 63 * return the answer in three pieces: 64 * y*log(x) = m*ln2 + hi + lo, 65 * where m is an integer. 66 * 3. Return x**y = exp(y*log(x)) 67 * = 2^m * ( exp(hi+lo) ). 68 * 69 * Special cases: 70 * (anything) ** 0 is 1 ; 71 * (anything) ** 1 is itself; 72 * (anything) ** NaN is NaN; 73 * NaN ** (anything except 0) is NaN; 74 * +(anything > 1) ** +INF is +INF; 75 * -(anything > 1) ** +INF is NaN; 76 * +-(anything > 1) ** -INF is +0; 77 * +-(anything < 1) ** +INF is +0; 78 * +(anything < 1) ** -INF is +INF; 79 * -(anything < 1) ** -INF is NaN; 80 * +-1 ** +-INF is NaN and signal INVALID; 81 * +0 ** +(anything except 0, NaN) is +0; 82 * -0 ** +(anything except 0, NaN, odd integer) is +0; 83 * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; 84 * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; 85 * -0 ** (odd integer) = -( +0 ** (odd integer) ); 86 * +INF ** +(anything except 0,NaN) is +INF; 87 * +INF ** -(anything except 0,NaN) is +0; 88 * -INF ** (odd integer) = -( +INF ** (odd integer) ); 89 * -INF ** (even integer) = ( +INF ** (even integer) ); 90 * -INF ** -(anything except integer,NaN) is NaN with signal; 91 * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); 92 * -(anything except 0) ** (non-integer) is NaN with signal; 93 * 94 * Accuracy: 95 * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, 96 * and a Zilog Z8000, 97 * pow(integer,integer) 98 * always returns the correct integer provided it is representable. 99 * In a test run with 100,000 random arguments with 0 < x, y < 20.0 100 * on a VAX, the maximum observed error was 1.79 ulps (units in the 101 * last place). 102 * 103 * Constants : 104 * The hexadecimal values are the intended ones for the following constants. 105 * The decimal values may be used, provided that the compiler will convert 106 * from decimal to binary accurately enough to produce the hexadecimal values 107 * shown. 108 */ 109 110 #include <errno.h> 111 #include <math.h> 112 113 #include "mathimpl.h" 114 115 #if (defined(__vax__) || defined(tahoe)) 116 #define TRUNC(x) x = (double) (float) x 117 #define _IEEE 0 118 #else 119 #define _IEEE 1 120 #define endian (((*(int *) &one)) ? 1 : 0) 121 #define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000 122 #define infnan(x) 0.0 123 #endif /* __vax__ or tahoe */ 124 125 static const double zero=0.0, one=1.0, two=2.0, negone= -1.0; 126 127 static double pow_P (double, double); 128 129 float 130 powf(float x, float y) 131 { 132 return pow((double) x, (double) (y)); 133 } 134 135 double 136 pow(double x, double y) 137 { 138 double t; 139 if (y==zero) 140 return (one); 141 else if (y==one || (_IEEE && x != x)) 142 return (x); /* if x is NaN or y=1 */ 143 else if (_IEEE && y!=y) /* if y is NaN */ 144 return (y); 145 else if (!finite(y)) /* if y is INF */ 146 if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */ 147 return (y - y); 148 else if (t>one) 149 return ((y<0)? zero : ((x<zero)? y-y : y)); 150 else 151 return ((y>0)? zero : ((x<0)? y-y : -y)); 152 else if (y==two) 153 return (x*x); 154 else if (y==negone) 155 return (one/x); 156 /* x > 0, x == +0 */ 157 else if (copysign(one, x) == one) 158 return (pow_P(x, y)); 159 160 /* sign(x)= -1 */ 161 /* if y is an even integer */ 162 else if ( (t=drem(y,two)) == zero) 163 return (pow_P(-x, y)); 164 165 /* if y is an odd integer */ 166 else if (copysign(t,one) == one) 167 return (-pow_P(-x, y)); 168 169 /* Henceforth y is not an integer */ 170 else if (x==zero) /* x is -0 */ 171 return ((y>zero)? -x : one/(-x)); 172 else if (_IEEE) 173 return (zero/zero); 174 else 175 return (infnan(EDOM)); 176 } 177 178 /* kernel function for x >= 0 */ 179 static double 180 pow_P(double x, double y) 181 { 182 struct Double s, t; 183 double huge = 1e300, tiny = 1e-300; 184 185 if (x == zero) { 186 if (y > zero) 187 return (zero); 188 else if (_IEEE) 189 return (huge*huge); 190 else 191 return (infnan(ERANGE)); 192 } 193 if (x == one) 194 return (one); 195 if (!finite(x)) { 196 if (y < zero) 197 return (zero); 198 else if (_IEEE) 199 return (huge*huge); 200 else 201 return (infnan(ERANGE)); 202 } 203 if (y >= 7e18) { /* infinity */ 204 if (x < 1) 205 return(tiny*tiny); 206 else if (_IEEE) 207 return (huge*huge); 208 else 209 return (infnan(ERANGE)); 210 } 211 212 /* Return exp(y*log(x)), using simulated extended */ 213 /* precision for the log and the multiply. */ 214 215 s = __log__D(x); 216 t.a = y; 217 TRUNC(t.a); 218 t.b = y - t.a; 219 t.b = s.b*y + t.b*s.a; 220 t.a *= s.a; 221 s.a = t.a + t.b; 222 s.b = (t.a - s.a) + t.b; 223 return (__exp__D(s.a, s.b)); 224 } 225