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