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