1 |
|
|
/* $OpenBSD: e_powl.c,v 1.7 2017/01/21 08:29:13 krw Exp $ */ |
2 |
|
|
|
3 |
|
|
/* |
4 |
|
|
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
5 |
|
|
* |
6 |
|
|
* Permission to use, copy, modify, and distribute this software for any |
7 |
|
|
* purpose with or without fee is hereby granted, provided that the above |
8 |
|
|
* copyright notice and this permission notice appear in all copies. |
9 |
|
|
* |
10 |
|
|
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
11 |
|
|
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
12 |
|
|
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR |
13 |
|
|
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
14 |
|
|
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN |
15 |
|
|
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF |
16 |
|
|
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. |
17 |
|
|
*/ |
18 |
|
|
|
19 |
|
|
/* powl.c |
20 |
|
|
* |
21 |
|
|
* Power function, long double precision |
22 |
|
|
* |
23 |
|
|
* |
24 |
|
|
* |
25 |
|
|
* SYNOPSIS: |
26 |
|
|
* |
27 |
|
|
* long double x, y, z, powl(); |
28 |
|
|
* |
29 |
|
|
* z = powl( x, y ); |
30 |
|
|
* |
31 |
|
|
* |
32 |
|
|
* |
33 |
|
|
* DESCRIPTION: |
34 |
|
|
* |
35 |
|
|
* Computes x raised to the yth power. Analytically, |
36 |
|
|
* |
37 |
|
|
* x**y = exp( y log(x) ). |
38 |
|
|
* |
39 |
|
|
* Following Cody and Waite, this program uses a lookup table |
40 |
|
|
* of 2**-i/32 and pseudo extended precision arithmetic to |
41 |
|
|
* obtain several extra bits of accuracy in both the logarithm |
42 |
|
|
* and the exponential. |
43 |
|
|
* |
44 |
|
|
* |
45 |
|
|
* |
46 |
|
|
* ACCURACY: |
47 |
|
|
* |
48 |
|
|
* The relative error of pow(x,y) can be estimated |
49 |
|
|
* by y dl ln(2), where dl is the absolute error of |
50 |
|
|
* the internally computed base 2 logarithm. At the ends |
51 |
|
|
* of the approximation interval the logarithm equal 1/32 |
52 |
|
|
* and its relative error is about 1 lsb = 1.1e-19. Hence |
53 |
|
|
* the predicted relative error in the result is 2.3e-21 y . |
54 |
|
|
* |
55 |
|
|
* Relative error: |
56 |
|
|
* arithmetic domain # trials peak rms |
57 |
|
|
* |
58 |
|
|
* IEEE +-1000 40000 2.8e-18 3.7e-19 |
59 |
|
|
* .001 < x < 1000, with log(x) uniformly distributed. |
60 |
|
|
* -1000 < y < 1000, y uniformly distributed. |
61 |
|
|
* |
62 |
|
|
* IEEE 0,8700 60000 6.5e-18 1.0e-18 |
63 |
|
|
* 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. |
64 |
|
|
* |
65 |
|
|
* |
66 |
|
|
* ERROR MESSAGES: |
67 |
|
|
* |
68 |
|
|
* message condition value returned |
69 |
|
|
* pow overflow x**y > MAXNUM INFINITY |
70 |
|
|
* pow underflow x**y < 1/MAXNUM 0.0 |
71 |
|
|
* pow domain x<0 and y noninteger 0.0 |
72 |
|
|
* |
73 |
|
|
*/ |
74 |
|
|
|
75 |
|
|
#include <float.h> |
76 |
|
|
#include <math.h> |
77 |
|
|
|
78 |
|
|
#include "math_private.h" |
79 |
|
|
|
80 |
|
|
/* Table size */ |
81 |
|
|
#define NXT 32 |
82 |
|
|
/* log2(Table size) */ |
83 |
|
|
#define LNXT 5 |
84 |
|
|
|
85 |
|
|
/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) |
86 |
|
|
* on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 |
87 |
|
|
*/ |
88 |
|
|
static long double P[] = { |
89 |
|
|
8.3319510773868690346226E-4L, |
90 |
|
|
4.9000050881978028599627E-1L, |
91 |
|
|
1.7500123722550302671919E0L, |
92 |
|
|
1.4000100839971580279335E0L, |
93 |
|
|
}; |
94 |
|
|
static long double Q[] = { |
95 |
|
|
/* 1.0000000000000000000000E0L,*/ |
96 |
|
|
5.2500282295834889175431E0L, |
97 |
|
|
8.4000598057587009834666E0L, |
98 |
|
|
4.2000302519914740834728E0L, |
99 |
|
|
}; |
100 |
|
|
/* A[i] = 2^(-i/32), rounded to IEEE long double precision. |
101 |
|
|
* If i is even, A[i] + B[i/2] gives additional accuracy. |
102 |
|
|
*/ |
103 |
|
|
static long double A[33] = { |
104 |
|
|
1.0000000000000000000000E0L, |
105 |
|
|
9.7857206208770013448287E-1L, |
106 |
|
|
9.5760328069857364691013E-1L, |
107 |
|
|
9.3708381705514995065011E-1L, |
108 |
|
|
9.1700404320467123175367E-1L, |
109 |
|
|
8.9735453750155359320742E-1L, |
110 |
|
|
8.7812608018664974155474E-1L, |
111 |
|
|
8.5930964906123895780165E-1L, |
112 |
|
|
8.4089641525371454301892E-1L, |
113 |
|
|
8.2287773907698242225554E-1L, |
114 |
|
|
8.0524516597462715409607E-1L, |
115 |
|
|
7.8799042255394324325455E-1L, |
116 |
|
|
7.7110541270397041179298E-1L, |
117 |
|
|
7.5458221379671136985669E-1L, |
118 |
|
|
7.3841307296974965571198E-1L, |
119 |
|
|
7.2259040348852331001267E-1L, |
120 |
|
|
7.0710678118654752438189E-1L, |
121 |
|
|
6.9195494098191597746178E-1L, |
122 |
|
|
6.7712777346844636413344E-1L, |
123 |
|
|
6.6261832157987064729696E-1L, |
124 |
|
|
6.4841977732550483296079E-1L, |
125 |
|
|
6.3452547859586661129850E-1L, |
126 |
|
|
6.2092890603674202431705E-1L, |
127 |
|
|
6.0762367999023443907803E-1L, |
128 |
|
|
5.9460355750136053334378E-1L, |
129 |
|
|
5.8186242938878875689693E-1L, |
130 |
|
|
5.6939431737834582684856E-1L, |
131 |
|
|
5.5719337129794626814472E-1L, |
132 |
|
|
5.4525386633262882960438E-1L, |
133 |
|
|
5.3357020033841180906486E-1L, |
134 |
|
|
5.2213689121370692017331E-1L, |
135 |
|
|
5.1094857432705833910408E-1L, |
136 |
|
|
5.0000000000000000000000E-1L, |
137 |
|
|
}; |
138 |
|
|
static long double B[17] = { |
139 |
|
|
0.0000000000000000000000E0L, |
140 |
|
|
2.6176170809902549338711E-20L, |
141 |
|
|
-1.0126791927256478897086E-20L, |
142 |
|
|
1.3438228172316276937655E-21L, |
143 |
|
|
1.2207982955417546912101E-20L, |
144 |
|
|
-6.3084814358060867200133E-21L, |
145 |
|
|
1.3164426894366316434230E-20L, |
146 |
|
|
-1.8527916071632873716786E-20L, |
147 |
|
|
1.8950325588932570796551E-20L, |
148 |
|
|
1.5564775779538780478155E-20L, |
149 |
|
|
6.0859793637556860974380E-21L, |
150 |
|
|
-2.0208749253662532228949E-20L, |
151 |
|
|
1.4966292219224761844552E-20L, |
152 |
|
|
3.3540909728056476875639E-21L, |
153 |
|
|
-8.6987564101742849540743E-22L, |
154 |
|
|
-1.2327176863327626135542E-20L, |
155 |
|
|
0.0000000000000000000000E0L, |
156 |
|
|
}; |
157 |
|
|
|
158 |
|
|
/* 2^x = 1 + x P(x), |
159 |
|
|
* on the interval -1/32 <= x <= 0 |
160 |
|
|
*/ |
161 |
|
|
static long double R[] = { |
162 |
|
|
1.5089970579127659901157E-5L, |
163 |
|
|
1.5402715328927013076125E-4L, |
164 |
|
|
1.3333556028915671091390E-3L, |
165 |
|
|
9.6181291046036762031786E-3L, |
166 |
|
|
5.5504108664798463044015E-2L, |
167 |
|
|
2.4022650695910062854352E-1L, |
168 |
|
|
6.9314718055994530931447E-1L, |
169 |
|
|
}; |
170 |
|
|
|
171 |
|
|
#define douba(k) A[k] |
172 |
|
|
#define doubb(k) B[k] |
173 |
|
|
#define MEXP (NXT*16384.0L) |
174 |
|
|
/* The following if denormal numbers are supported, else -MEXP: */ |
175 |
|
|
#define MNEXP (-NXT*(16384.0L+64.0L)) |
176 |
|
|
/* log2(e) - 1 */ |
177 |
|
|
#define LOG2EA 0.44269504088896340735992L |
178 |
|
|
|
179 |
|
|
#define F W |
180 |
|
|
#define Fa Wa |
181 |
|
|
#define Fb Wb |
182 |
|
|
#define G W |
183 |
|
|
#define Ga Wa |
184 |
|
|
#define Gb u |
185 |
|
|
#define H W |
186 |
|
|
#define Ha Wb |
187 |
|
|
#define Hb Wb |
188 |
|
|
|
189 |
|
|
static const long double MAXLOGL = 1.1356523406294143949492E4L; |
190 |
|
|
static const long double MINLOGL = -1.13994985314888605586758E4L; |
191 |
|
|
static const long double LOGE2L = 6.9314718055994530941723E-1L; |
192 |
|
|
static volatile long double z; |
193 |
|
|
static long double w, W, Wa, Wb, ya, yb, u; |
194 |
|
|
static const long double huge = 0x1p10000L; |
195 |
|
|
#if 0 /* XXX Prevent gcc from erroneously constant folding this. */ |
196 |
|
|
static const long double twom10000 = 0x1p-10000L; |
197 |
|
|
#else |
198 |
|
|
static volatile long double twom10000 = 0x1p-10000L; |
199 |
|
|
#endif |
200 |
|
|
|
201 |
|
|
static long double reducl( long double ); |
202 |
|
|
static long double powil ( long double, int ); |
203 |
|
|
|
204 |
|
|
long double |
205 |
|
|
powl(long double x, long double y) |
206 |
|
|
{ |
207 |
|
|
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */ |
208 |
|
|
int i, nflg, iyflg, yoddint; |
209 |
|
|
long e; |
210 |
|
|
|
211 |
|
|
if( y == 0.0L ) |
212 |
|
|
return( 1.0L ); |
213 |
|
|
|
214 |
|
|
if( x == 1.0L ) |
215 |
|
|
return( 1.0L ); |
216 |
|
|
|
217 |
|
|
if( isnan(x) ) |
218 |
|
|
return( x ); |
219 |
|
|
if( isnan(y) ) |
220 |
|
|
return( y ); |
221 |
|
|
|
222 |
|
|
if( y == 1.0L ) |
223 |
|
|
return( x ); |
224 |
|
|
|
225 |
|
|
if( !isfinite(y) && x == -1.0L ) |
226 |
|
|
return( 1.0L ); |
227 |
|
|
|
228 |
|
|
if( y >= LDBL_MAX ) |
229 |
|
|
{ |
230 |
|
|
if( x > 1.0L ) |
231 |
|
|
return( INFINITY ); |
232 |
|
|
if( x > 0.0L && x < 1.0L ) |
233 |
|
|
return( 0.0L ); |
234 |
|
|
if( x < -1.0L ) |
235 |
|
|
return( INFINITY ); |
236 |
|
|
if( x > -1.0L && x < 0.0L ) |
237 |
|
|
return( 0.0L ); |
238 |
|
|
} |
239 |
|
|
if( y <= -LDBL_MAX ) |
240 |
|
|
{ |
241 |
|
|
if( x > 1.0L ) |
242 |
|
|
return( 0.0L ); |
243 |
|
|
if( x > 0.0L && x < 1.0L ) |
244 |
|
|
return( INFINITY ); |
245 |
|
|
if( x < -1.0L ) |
246 |
|
|
return( 0.0L ); |
247 |
|
|
if( x > -1.0L && x < 0.0L ) |
248 |
|
|
return( INFINITY ); |
249 |
|
|
} |
250 |
|
|
if( x >= LDBL_MAX ) |
251 |
|
|
{ |
252 |
|
|
if( y > 0.0L ) |
253 |
|
|
return( INFINITY ); |
254 |
|
|
return( 0.0L ); |
255 |
|
|
} |
256 |
|
|
|
257 |
|
|
w = floorl(y); |
258 |
|
|
/* Set iyflg to 1 if y is an integer. */ |
259 |
|
|
iyflg = 0; |
260 |
|
|
if( w == y ) |
261 |
|
|
iyflg = 1; |
262 |
|
|
|
263 |
|
|
/* Test for odd integer y. */ |
264 |
|
|
yoddint = 0; |
265 |
|
|
if( iyflg ) |
266 |
|
|
{ |
267 |
|
|
ya = fabsl(y); |
268 |
|
|
ya = floorl(0.5L * ya); |
269 |
|
|
yb = 0.5L * fabsl(w); |
270 |
|
|
if( ya != yb ) |
271 |
|
|
yoddint = 1; |
272 |
|
|
} |
273 |
|
|
|
274 |
|
|
if( x <= -LDBL_MAX ) |
275 |
|
|
{ |
276 |
|
|
if( y > 0.0L ) |
277 |
|
|
{ |
278 |
|
|
if( yoddint ) |
279 |
|
|
return( -INFINITY ); |
280 |
|
|
return( INFINITY ); |
281 |
|
|
} |
282 |
|
|
if( y < 0.0L ) |
283 |
|
|
{ |
284 |
|
|
if( yoddint ) |
285 |
|
|
return( -0.0L ); |
286 |
|
|
return( 0.0 ); |
287 |
|
|
} |
288 |
|
|
} |
289 |
|
|
|
290 |
|
|
|
291 |
|
|
nflg = 0; /* flag = 1 if x<0 raised to integer power */ |
292 |
|
|
if( x <= 0.0L ) |
293 |
|
|
{ |
294 |
|
|
if( x == 0.0L ) |
295 |
|
|
{ |
296 |
|
|
if( y < 0.0 ) |
297 |
|
|
{ |
298 |
|
|
if( signbit(x) && yoddint ) |
299 |
|
|
return( -INFINITY ); |
300 |
|
|
return( INFINITY ); |
301 |
|
|
} |
302 |
|
|
if( y > 0.0 ) |
303 |
|
|
{ |
304 |
|
|
if( signbit(x) && yoddint ) |
305 |
|
|
return( -0.0L ); |
306 |
|
|
return( 0.0 ); |
307 |
|
|
} |
308 |
|
|
if( y == 0.0L ) |
309 |
|
|
return( 1.0L ); /* 0**0 */ |
310 |
|
|
else |
311 |
|
|
return( 0.0L ); /* 0**y */ |
312 |
|
|
} |
313 |
|
|
else |
314 |
|
|
{ |
315 |
|
|
if( iyflg == 0 ) |
316 |
|
|
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ |
317 |
|
|
nflg = 1; |
318 |
|
|
} |
319 |
|
|
} |
320 |
|
|
|
321 |
|
|
/* Integer power of an integer. */ |
322 |
|
|
|
323 |
|
|
if( iyflg ) |
324 |
|
|
{ |
325 |
|
|
i = w; |
326 |
|
|
w = floorl(x); |
327 |
|
|
if( (w == x) && (fabsl(y) < 32768.0) ) |
328 |
|
|
{ |
329 |
|
|
w = powil( x, (int) y ); |
330 |
|
|
return( w ); |
331 |
|
|
} |
332 |
|
|
} |
333 |
|
|
|
334 |
|
|
|
335 |
|
|
if( nflg ) |
336 |
|
|
x = fabsl(x); |
337 |
|
|
|
338 |
|
|
/* separate significand from exponent */ |
339 |
|
|
x = frexpl( x, &i ); |
340 |
|
|
e = i; |
341 |
|
|
|
342 |
|
|
/* find significand in antilog table A[] */ |
343 |
|
|
i = 1; |
344 |
|
|
if( x <= douba(17) ) |
345 |
|
|
i = 17; |
346 |
|
|
if( x <= douba(i+8) ) |
347 |
|
|
i += 8; |
348 |
|
|
if( x <= douba(i+4) ) |
349 |
|
|
i += 4; |
350 |
|
|
if( x <= douba(i+2) ) |
351 |
|
|
i += 2; |
352 |
|
|
if( x >= douba(1) ) |
353 |
|
|
i = -1; |
354 |
|
|
i += 1; |
355 |
|
|
|
356 |
|
|
|
357 |
|
|
/* Find (x - A[i])/A[i] |
358 |
|
|
* in order to compute log(x/A[i]): |
359 |
|
|
* |
360 |
|
|
* log(x) = log( a x/a ) = log(a) + log(x/a) |
361 |
|
|
* |
362 |
|
|
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a |
363 |
|
|
*/ |
364 |
|
|
x -= douba(i); |
365 |
|
|
x -= doubb(i/2); |
366 |
|
|
x /= douba(i); |
367 |
|
|
|
368 |
|
|
|
369 |
|
|
/* rational approximation for log(1+v): |
370 |
|
|
* |
371 |
|
|
* log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) |
372 |
|
|
*/ |
373 |
|
|
z = x*x; |
374 |
|
|
w = x * ( z * __polevll( x, P, 3 ) / __p1evll( x, Q, 3 ) ); |
375 |
|
|
w = w - ldexpl( z, -1 ); /* w - 0.5 * z */ |
376 |
|
|
|
377 |
|
|
/* Convert to base 2 logarithm: |
378 |
|
|
* multiply by log2(e) = 1 + LOG2EA |
379 |
|
|
*/ |
380 |
|
|
z = LOG2EA * w; |
381 |
|
|
z += w; |
382 |
|
|
z += LOG2EA * x; |
383 |
|
|
z += x; |
384 |
|
|
|
385 |
|
|
/* Compute exponent term of the base 2 logarithm. */ |
386 |
|
|
w = -i; |
387 |
|
|
w = ldexpl( w, -LNXT ); /* divide by NXT */ |
388 |
|
|
w += e; |
389 |
|
|
/* Now base 2 log of x is w + z. */ |
390 |
|
|
|
391 |
|
|
/* Multiply base 2 log by y, in extended precision. */ |
392 |
|
|
|
393 |
|
|
/* separate y into large part ya |
394 |
|
|
* and small part yb less than 1/NXT |
395 |
|
|
*/ |
396 |
|
|
ya = reducl(y); |
397 |
|
|
yb = y - ya; |
398 |
|
|
|
399 |
|
|
/* (w+z)(ya+yb) |
400 |
|
|
* = w*ya + w*yb + z*y |
401 |
|
|
*/ |
402 |
|
|
F = z * y + w * yb; |
403 |
|
|
Fa = reducl(F); |
404 |
|
|
Fb = F - Fa; |
405 |
|
|
|
406 |
|
|
G = Fa + w * ya; |
407 |
|
|
Ga = reducl(G); |
408 |
|
|
Gb = G - Ga; |
409 |
|
|
|
410 |
|
|
H = Fb + Gb; |
411 |
|
|
Ha = reducl(H); |
412 |
|
|
w = ldexpl( Ga+Ha, LNXT ); |
413 |
|
|
|
414 |
|
|
/* Test the power of 2 for overflow */ |
415 |
|
|
if( w > MEXP ) |
416 |
|
|
return (huge * huge); /* overflow */ |
417 |
|
|
|
418 |
|
|
if( w < MNEXP ) |
419 |
|
|
return (twom10000 * twom10000); /* underflow */ |
420 |
|
|
|
421 |
|
|
e = w; |
422 |
|
|
Hb = H - Ha; |
423 |
|
|
|
424 |
|
|
if( Hb > 0.0L ) |
425 |
|
|
{ |
426 |
|
|
e += 1; |
427 |
|
|
Hb -= (1.0L/NXT); /*0.0625L;*/ |
428 |
|
|
} |
429 |
|
|
|
430 |
|
|
/* Now the product y * log2(x) = Hb + e/NXT. |
431 |
|
|
* |
432 |
|
|
* Compute base 2 exponential of Hb, |
433 |
|
|
* where -0.0625 <= Hb <= 0. |
434 |
|
|
*/ |
435 |
|
|
z = Hb * __polevll( Hb, R, 6 ); /* z = 2**Hb - 1 */ |
436 |
|
|
|
437 |
|
|
/* Express e/NXT as an integer plus a negative number of (1/NXT)ths. |
438 |
|
|
* Find lookup table entry for the fractional power of 2. |
439 |
|
|
*/ |
440 |
|
|
if( e < 0 ) |
441 |
|
|
i = 0; |
442 |
|
|
else |
443 |
|
|
i = 1; |
444 |
|
|
i = e/NXT + i; |
445 |
|
|
e = NXT*i - e; |
446 |
|
|
w = douba( e ); |
447 |
|
|
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ |
448 |
|
|
z = z + w; |
449 |
|
|
z = ldexpl( z, i ); /* multiply by integer power of 2 */ |
450 |
|
|
|
451 |
|
|
if( nflg ) |
452 |
|
|
{ |
453 |
|
|
/* For negative x, |
454 |
|
|
* find out if the integer exponent |
455 |
|
|
* is odd or even. |
456 |
|
|
*/ |
457 |
|
|
w = ldexpl( y, -1 ); |
458 |
|
|
w = floorl(w); |
459 |
|
|
w = ldexpl( w, 1 ); |
460 |
|
|
if( w != y ) |
461 |
|
|
z = -z; /* odd exponent */ |
462 |
|
|
} |
463 |
|
|
|
464 |
|
|
return( z ); |
465 |
|
|
} |
466 |
|
|
DEF_STD(powl); |
467 |
|
|
|
468 |
|
|
|
469 |
|
|
/* Find a multiple of 1/NXT that is within 1/NXT of x. */ |
470 |
|
|
static long double |
471 |
|
|
reducl(long double x) |
472 |
|
|
{ |
473 |
|
|
long double t; |
474 |
|
|
|
475 |
|
|
t = ldexpl( x, LNXT ); |
476 |
|
|
t = floorl( t ); |
477 |
|
|
t = ldexpl( t, -LNXT ); |
478 |
|
|
return(t); |
479 |
|
|
} |
480 |
|
|
|
481 |
|
|
/* powil.c |
482 |
|
|
* |
483 |
|
|
* Real raised to integer power, long double precision |
484 |
|
|
* |
485 |
|
|
* |
486 |
|
|
* |
487 |
|
|
* SYNOPSIS: |
488 |
|
|
* |
489 |
|
|
* long double x, y, powil(); |
490 |
|
|
* int n; |
491 |
|
|
* |
492 |
|
|
* y = powil( x, n ); |
493 |
|
|
* |
494 |
|
|
* |
495 |
|
|
* |
496 |
|
|
* DESCRIPTION: |
497 |
|
|
* |
498 |
|
|
* Returns argument x raised to the nth power. |
499 |
|
|
* The routine efficiently decomposes n as a sum of powers of |
500 |
|
|
* two. The desired power is a product of two-to-the-kth |
501 |
|
|
* powers of x. Thus to compute the 32767 power of x requires |
502 |
|
|
* 28 multiplications instead of 32767 multiplications. |
503 |
|
|
* |
504 |
|
|
* |
505 |
|
|
* |
506 |
|
|
* ACCURACY: |
507 |
|
|
* |
508 |
|
|
* |
509 |
|
|
* Relative error: |
510 |
|
|
* arithmetic x domain n domain # trials peak rms |
511 |
|
|
* IEEE .001,1000 -1022,1023 50000 4.3e-17 7.8e-18 |
512 |
|
|
* IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18 |
513 |
|
|
* IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17 |
514 |
|
|
* |
515 |
|
|
* Returns MAXNUM on overflow, zero on underflow. |
516 |
|
|
* |
517 |
|
|
*/ |
518 |
|
|
|
519 |
|
|
static long double |
520 |
|
|
powil(long double x, int nn) |
521 |
|
|
{ |
522 |
|
|
long double ww, y; |
523 |
|
|
long double s; |
524 |
|
|
int n, e, sign, asign, lx; |
525 |
|
|
|
526 |
|
|
if( x == 0.0L ) |
527 |
|
|
{ |
528 |
|
|
if( nn == 0 ) |
529 |
|
|
return( 1.0L ); |
530 |
|
|
else if( nn < 0 ) |
531 |
|
|
return( LDBL_MAX ); |
532 |
|
|
else |
533 |
|
|
return( 0.0L ); |
534 |
|
|
} |
535 |
|
|
|
536 |
|
|
if( nn == 0 ) |
537 |
|
|
return( 1.0L ); |
538 |
|
|
|
539 |
|
|
|
540 |
|
|
if( x < 0.0L ) |
541 |
|
|
{ |
542 |
|
|
asign = -1; |
543 |
|
|
x = -x; |
544 |
|
|
} |
545 |
|
|
else |
546 |
|
|
asign = 0; |
547 |
|
|
|
548 |
|
|
|
549 |
|
|
if( nn < 0 ) |
550 |
|
|
{ |
551 |
|
|
sign = -1; |
552 |
|
|
n = -nn; |
553 |
|
|
} |
554 |
|
|
else |
555 |
|
|
{ |
556 |
|
|
sign = 1; |
557 |
|
|
n = nn; |
558 |
|
|
} |
559 |
|
|
|
560 |
|
|
/* Overflow detection */ |
561 |
|
|
|
562 |
|
|
/* Calculate approximate logarithm of answer */ |
563 |
|
|
s = x; |
564 |
|
|
s = frexpl( s, &lx ); |
565 |
|
|
e = (lx - 1)*n; |
566 |
|
|
if( (e == 0) || (e > 64) || (e < -64) ) |
567 |
|
|
{ |
568 |
|
|
s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L); |
569 |
|
|
s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L; |
570 |
|
|
} |
571 |
|
|
else |
572 |
|
|
{ |
573 |
|
|
s = LOGE2L * e; |
574 |
|
|
} |
575 |
|
|
|
576 |
|
|
if( s > MAXLOGL ) |
577 |
|
|
return (huge * huge); /* overflow */ |
578 |
|
|
|
579 |
|
|
if( s < MINLOGL ) |
580 |
|
|
return (twom10000 * twom10000); /* underflow */ |
581 |
|
|
/* Handle tiny denormal answer, but with less accuracy |
582 |
|
|
* since roundoff error in 1.0/x will be amplified. |
583 |
|
|
* The precise demarcation should be the gradual underflow threshold. |
584 |
|
|
*/ |
585 |
|
|
if( s < (-MAXLOGL+2.0L) ) |
586 |
|
|
{ |
587 |
|
|
x = 1.0L/x; |
588 |
|
|
sign = -sign; |
589 |
|
|
} |
590 |
|
|
|
591 |
|
|
/* First bit of the power */ |
592 |
|
|
if( n & 1 ) |
593 |
|
|
y = x; |
594 |
|
|
|
595 |
|
|
else |
596 |
|
|
{ |
597 |
|
|
y = 1.0L; |
598 |
|
|
asign = 0; |
599 |
|
|
} |
600 |
|
|
|
601 |
|
|
ww = x; |
602 |
|
|
n >>= 1; |
603 |
|
|
while( n ) |
604 |
|
|
{ |
605 |
|
|
ww = ww * ww; /* arg to the 2-to-the-kth power */ |
606 |
|
|
if( n & 1 ) /* if that bit is set, then include in product */ |
607 |
|
|
y *= ww; |
608 |
|
|
n >>= 1; |
609 |
|
|
} |
610 |
|
|
|
611 |
|
|
if( asign ) |
612 |
|
|
y = -y; /* odd power of negative number */ |
613 |
|
|
if( sign < 0 ) |
614 |
|
|
y = 1.0L/y; |
615 |
|
|
return(y); |
616 |
|
|
} |