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
pow(x,y)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
pow_P(double x,double y)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