SDL 2.0
e_pow.c
Go to the documentation of this file.
1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/* __ieee754_pow(x,y) return x**y
13 *
14 * n
15 * Method: Let x = 2 * (1+f)
16 * 1. Compute and return log2(x) in two pieces:
17 * log2(x) = w1 + w2,
18 * where w1 has 53-24 = 29 bit trailing zeros.
19 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
20 * arithmetic, where |y'|<=0.5.
21 * 3. Return x**y = 2**n*exp(y'*log2)
22 *
23 * Special cases:
24 * 1. +-1 ** anything is 1.0
25 * 2. +-1 ** +-INF is 1.0
26 * 3. (anything) ** 0 is 1
27 * 4. (anything) ** 1 is itself
28 * 5. (anything) ** NAN is NAN
29 * 6. NAN ** (anything except 0) is NAN
30 * 7. +-(|x| > 1) ** +INF is +INF
31 * 8. +-(|x| > 1) ** -INF is +0
32 * 9. +-(|x| < 1) ** +INF is +0
33 * 10 +-(|x| < 1) ** -INF is +INF
34 * 11. +0 ** (+anything except 0, NAN) is +0
35 * 12. -0 ** (+anything except 0, NAN, odd integer) is +0
36 * 13. +0 ** (-anything except 0, NAN) is +INF
37 * 14. -0 ** (-anything except 0, NAN, odd integer) is +INF
38 * 15. -0 ** (odd integer) = -( +0 ** (odd integer) )
39 * 16. +INF ** (+anything except 0,NAN) is +INF
40 * 17. +INF ** (-anything except 0,NAN) is +0
41 * 18. -INF ** (anything) = -0 ** (-anything)
42 * 19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
43 * 20. (-anything except 0 and inf) ** (non-integer) is NAN
44 *
45 * Accuracy:
46 * pow(x,y) returns x**y nearly rounded. In particular
47 * pow(integer,integer)
48 * always returns the correct integer provided it is
49 * representable.
50 *
51 * Constants :
52 * The hexadecimal values are the intended ones for the following
53 * constants. The decimal values may be used, provided that the
54 * compiler will convert from decimal to binary accurately enough
55 * to produce the hexadecimal values shown.
56 */
57
58#include "math_libm.h"
59#include "math_private.h"
60
61#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
62/* C4756: overflow in constant arithmetic */
63#pragma warning ( disable : 4756 )
64#endif
65
66#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
67#undef huge
68#endif
69
70static const double
71bp[] = {1.0, 1.5,},
72dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
73dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
74zero = 0.0,
75one = 1.0,
76two = 2.0,
77two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
78huge = 1.0e300,
79tiny = 1.0e-300,
80 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
81L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
82L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
83L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
84L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
85L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
86L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
87P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
88P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
89P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
90P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
91P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
92lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
93lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
94lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
95ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
96cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
97cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
98cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
99ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
100ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
101ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
102
103double attribute_hidden __ieee754_pow(double x, double y)
104{
105 double z,ax,z_h,z_l,p_h,p_l;
106 double y1,t1,t2,r,s,t,u,v,w;
107 int32_t i,j,k,yisint,n;
108 int32_t hx,hy,ix,iy;
109 u_int32_t lx,ly;
110
111 EXTRACT_WORDS(hx,lx,x);
112 /* x==1: 1**y = 1 (even if y is NaN) */
113 if (hx==0x3ff00000 && lx==0) {
114 return x;
115 }
116 ix = hx&0x7fffffff;
117
118 EXTRACT_WORDS(hy,ly,y);
119 iy = hy&0x7fffffff;
120
121 /* y==zero: x**0 = 1 */
122 if((iy|ly)==0) return one;
123
124 /* +-NaN return x+y */
125 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
126 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
127 return x+y;
128
129 /* determine if y is an odd int when x < 0
130 * yisint = 0 ... y is not an integer
131 * yisint = 1 ... y is an odd int
132 * yisint = 2 ... y is an even int
133 */
134 yisint = 0;
135 if(hx<0) {
136 if(iy>=0x43400000) yisint = 2; /* even integer y */
137 else if(iy>=0x3ff00000) {
138 k = (iy>>20)-0x3ff; /* exponent */
139 if(k>20) {
140 j = ly>>(52-k);
141 if((j<<(52-k))==ly) yisint = 2-(j&1);
142 } else if(ly==0) {
143 j = iy>>(20-k);
144 if((j<<(20-k))==iy) yisint = 2-(j&1);
145 }
146 }
147 }
148
149 /* special value of y */
150 if(ly==0) {
151 if (iy==0x7ff00000) { /* y is +-inf */
152 if (((ix-0x3ff00000)|lx)==0)
153 return one; /* +-1**+-inf is 1 (yes, weird rule) */
154 if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
155 return (hy>=0) ? y : zero;
156 /* (|x|<1)**-,+inf = inf,0 */
157 return (hy<0) ? -y : zero;
158 }
159 if(iy==0x3ff00000) { /* y is +-1 */
160 if(hy<0) return one/x; else return x;
161 }
162 if(hy==0x40000000) return x*x; /* y is 2 */
163 if(hy==0x3fe00000) { /* y is 0.5 */
164 if(hx>=0) /* x >= +0 */
165 return __ieee754_sqrt(x);
166 }
167 }
168
169 ax = fabs(x);
170 /* special value of x */
171 if(lx==0) {
172 if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
173 z = ax; /*x is +-0,+-inf,+-1*/
174 if(hy<0) z = one/z; /* z = (1/|x|) */
175 if(hx<0) {
176 if(((ix-0x3ff00000)|yisint)==0) {
177 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
178 } else if(yisint==1)
179 z = -z; /* (x<0)**odd = -(|x|**odd) */
180 }
181 return z;
182 }
183 }
184
185 /* (x<0)**(non-int) is NaN */
186 if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
187
188 /* |y| is huge */
189 if(iy>0x41e00000) { /* if |y| > 2**31 */
190 if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
191 if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
192 if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
193 }
194 /* over/underflow if x is not close to one */
195 if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
196 if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
197 /* now |1-x| is tiny <= 2**-20, suffice to compute
198 log(x) by x-x^2/2+x^3/3-x^4/4 */
199 t = x-1; /* t has 20 trailing zeros */
200 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
201 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
202 v = t*ivln2_l-w*ivln2;
203 t1 = u+v;
204 SET_LOW_WORD(t1,0);
205 t2 = v-(t1-u);
206 } else {
207 double s2,s_h,s_l,t_h,t_l;
208 n = 0;
209 /* take care subnormal number */
210 if(ix<0x00100000)
211 {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
212 n += ((ix)>>20)-0x3ff;
213 j = ix&0x000fffff;
214 /* determine interval */
215 ix = j|0x3ff00000; /* normalize ix */
216 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
217 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
218 else {k=0;n+=1;ix -= 0x00100000;}
219 SET_HIGH_WORD(ax,ix);
220
221 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
222 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
223 v = one/(ax+bp[k]);
224 s = u*v;
225 s_h = s;
226 SET_LOW_WORD(s_h,0);
227 /* t_h=ax+bp[k] High */
228 t_h = zero;
229 SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
230 t_l = ax - (t_h-bp[k]);
231 s_l = v*((u-s_h*t_h)-s_h*t_l);
232 /* compute log(ax) */
233 s2 = s*s;
234 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
235 r += s_l*(s_h+s);
236 s2 = s_h*s_h;
237 t_h = 3.0+s2+r;
238 SET_LOW_WORD(t_h,0);
239 t_l = r-((t_h-3.0)-s2);
240 /* u+v = s*(1+...) */
241 u = s_h*t_h;
242 v = s_l*t_h+t_l*s;
243 /* 2/(3log2)*(s+...) */
244 p_h = u+v;
245 SET_LOW_WORD(p_h,0);
246 p_l = v-(p_h-u);
247 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
248 z_l = cp_l*p_h+p_l*cp+dp_l[k];
249 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
250 t = (double)n;
251 t1 = (((z_h+z_l)+dp_h[k])+t);
252 SET_LOW_WORD(t1,0);
253 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
254 }
255
256 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
257 if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
258 s = -one;/* (-ve)**(odd int) */
259
260 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
261 y1 = y;
262 SET_LOW_WORD(y1,0);
263 p_l = (y-y1)*t1+y*t2;
264 p_h = y1*t1;
265 z = p_l+p_h;
267 if (j>=0x40900000) { /* z >= 1024 */
268 if(((j-0x40900000)|i)!=0) /* if z > 1024 */
269 return s*huge*huge; /* overflow */
270 else {
271 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
272 }
273 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
274 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
275 return s*tiny*tiny; /* underflow */
276 else {
277 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
278 }
279 }
280 /*
281 * compute 2**(p_h+p_l)
282 */
283 i = j&0x7fffffff;
284 k = (i>>20)-0x3ff;
285 n = 0;
286 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
287 n = j+(0x00100000>>(k+1));
288 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
289 t = zero;
290 SET_HIGH_WORD(t,n&~(0x000fffff>>k));
291 n = ((n&0x000fffff)|0x00100000)>>(20-k);
292 if(j<0) n = -n;
293 p_h -= t;
294 }
295 t = p_l+p_h;
296 SET_LOW_WORD(t,0);
297 u = t*lg2_h;
298 v = (p_l-(t-p_h))*lg2+t*lg2_l;
299 z = u+v;
300 w = v-(z-u);
301 t = z*z;
302 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
303 r = (z*t1)/(t1-two)-(w+z*w);
304 z = one-(r-z);
306 j += (n<<20);
307 if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
308 else SET_HIGH_WORD(z,j);
309 return s*z;
310}
311
312/*
313 * wrapper pow(x,y) return x**y
314 */
315#ifndef _IEEE_LIBM
316double pow(double x, double y)
317{
318 double z = __ieee754_pow(x, y);
319 if (_LIB_VERSION == _IEEE_|| isnan(y))
320 return z;
321 if (isnan(x)) {
322 if (y == 0.0)
323 return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */
324 return z;
325 }
326 if (x == 0.0) {
327 if (y == 0.0)
328 return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */
329 if (isfinite(y) && y < 0.0)
330 return __kernel_standard(x,y,23); /* pow(0.0,negative) */
331 return z;
332 }
333 if (!isfinite(z)) {
334 if (isfinite(x) && isfinite(y)) {
335 if (isnan(z))
336 return __kernel_standard(x, y, 24); /* pow neg**non-int */
337 return __kernel_standard(x, y, 21); /* pow overflow */
338 }
339 }
340 if (z == 0.0 && isfinite(x) && isfinite(y))
341 return __kernel_standard(x, y, 22); /* pow underflow */
342 return z;
343}
344#else
346#endif
signed int int32_t
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
const GLdouble * v
Definition: SDL_opengl.h:2064
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
GLdouble s
Definition: SDL_opengl.h:2063
GLfixed y1
GLdouble GLdouble z
GLdouble n
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
GLubyte GLubyte GLubyte GLubyte w
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int int in j)
Definition: SDL_x11sym.h:50
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int in i)
Definition: SDL_x11sym.h:50
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int int int return Display Window Cursor return Display Window return Display Drawable GC int int unsigned int unsigned int return Display Drawable GC int int _Xconst char int return Display Drawable GC int int unsigned int unsigned int return Display return Display Cursor return Display GC return XModifierKeymap return char Display Window int return Display return Display Atom return Display Window XWindowAttributes return Display Window return Display XEvent Bool(*) XPointer return Display Window Bool unsigned int int int Window Cursor Time return Display Window int return KeySym return Display _Xconst char Bool return Display _Xconst char return XKeyEvent char int KeySym XComposeStatus return Display int int int XVisualInfo return Display Window int int return _Xconst char return Display XEvent return Display Drawable GC XImage int int int int unsigned int unsigned int return Display Window Window Window int int int int unsigned int return Display Window Window int int return Display Window unsigned int unsigned int return Display Window Bool long XEvent return Display GC unsigned long return Display Window int Time return Display Window Window return Display Window unsigned long return Display Window XSizeHints Display Colormap XColor int return char int XTextProperty return XFontStruct _Xconst char int int int int XCharStruct return Display Window return Display Time return Display Colormap return Display Window Window int int unsigned int unsigned int int int return Display Window int return XExtensionInfo Display char XExtensionHooks int XPointer return XExtensionInfo XExtensionInfo Display return Display return Display unsigned long Display GC Display char long Display xReply int Bool return Display Bool return Display int SDL_X11_XESetEventToWireRetType return Display Window Window Window Window unsigned int return Display XShmSegmentInfo return Display Drawable GC XImage int int int int unsigned int unsigned int Boo k)
Definition: SDL_x11sym.h:213
double attribute_hidden __ieee754_pow(double x, double y)
Definition: e_pow.c:103
static const double lg2_h
Definition: e_pow.c:93
static const double cp_h
Definition: e_pow.c:97
static const double L2
Definition: e_pow.c:82
static const double L3
Definition: e_pow.c:83
static const double ovt
Definition: e_pow.c:95
static const double cp
Definition: e_pow.c:96
static const double ivln2
Definition: e_pow.c:99
static const double two
Definition: e_pow.c:76
static const double dp_h[]
Definition: e_pow.c:72
static const double tiny
Definition: e_pow.c:79
static const double zero
Definition: e_pow.c:74
static const double P1
Definition: e_pow.c:87
static const double L4
Definition: e_pow.c:84
static const double ivln2_h
Definition: e_pow.c:100
static const double two53
Definition: e_pow.c:77
static const double L1
Definition: e_pow.c:81
static const double bp[]
Definition: e_pow.c:71
static const double L6
Definition: e_pow.c:86
static const double L5
Definition: e_pow.c:85
static const double P3
Definition: e_pow.c:89
static const double cp_l
Definition: e_pow.c:98
static const double lg2_l
Definition: e_pow.c:94
static const double ivln2_l
Definition: e_pow.c:101
static const double huge
Definition: e_pow.c:78
static const double P2
Definition: e_pow.c:88
static const double dp_l[]
Definition: e_pow.c:73
static const double one
Definition: e_pow.c:75
static const double P4
Definition: e_pow.c:90
static const double lg2
Definition: e_pow.c:92
static const double P5
Definition: e_pow.c:91
unsigned int u_int32_t
Definition: math_private.h:31
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:99
#define scalbn
Definition: math_private.h:46
#define __ieee754_sqrt
Definition: math_private.h:48
#define strong_alias(x, y)
Definition: math_private.h:28
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
#define SET_LOW_WORD(d, v)
Definition: math_private.h:147
#define attribute_hidden
Definition: math_private.h:25
double fabs(double x)
Definition: s_fabs.c:22
libm_hidden_def(scalbln)
Definition: s_scalbn.c:66