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 
70 static const double
71 bp[] = {1.0, 1.5,},
72 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
73 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
74 zero = 0.0,
75 one = 1.0,
76 two = 2.0,
77 two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
78 huge = 1.0e300,
79 tiny = 1.0e-300,
80  /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
81 L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
82 L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
83 L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
84 L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
85 L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
86 L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
87 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
88 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
89 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
90 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
91 P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
92 lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
93 lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
94 lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
95 ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
96 cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
97 cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
98 cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
99 ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
100 ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
101 ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
102 
103 double 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;
266  EXTRACT_WORDS(j,i,z);
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);
305  GET_HIGH_WORD(j,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
316 double 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
347 libm_hidden_def(pow)
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
double attribute_hidden __ieee754_pow(double x, double y)
Definition: e_pow.c:103
static const double two53
Definition: e_pow.c:77
GLdouble GLdouble z
static const double dp_l[]
Definition: e_pow.c:73
GLdouble s
Definition: SDL_opengl.h:2063
const GLdouble * v
Definition: SDL_opengl.h:2064
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
static const double ivln2_h
Definition: e_pow.c:100
static const double P3
Definition: e_pow.c:89
static const double L1
Definition: e_pow.c:81
static const double L6
Definition: e_pow.c:86
static const double lg2
Definition: e_pow.c:92
static const double lg2_l
Definition: e_pow.c:94
#define strong_alias(x, y)
Definition: math_private.h:28
static const double ovt
Definition: e_pow.c:95
static const double tiny
Definition: e_pow.c:79
static const double bp[]
Definition: e_pow.c:71
#define attribute_hidden
Definition: math_private.h:25
#define __ieee754_sqrt
Definition: math_private.h:48
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
#define scalbn
Definition: math_private.h:46
static const double dp_h[]
Definition: e_pow.c:72
GLfixed y1
static const double L3
Definition: e_pow.c:83
static const double P1
Definition: e_pow.c:87
static const double P4
Definition: e_pow.c:90
unsigned int u_int32_t
Definition: math_private.h:31
static const double L5
Definition: e_pow.c:85
static const double two
Definition: e_pow.c:76
#define SET_LOW_WORD(d, v)
Definition: math_private.h:147
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:99
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
static const double one
Definition: e_pow.c:75
GLubyte GLubyte GLubyte GLubyte w
static const double lg2_h
Definition: e_pow.c:93
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
static const double cp_l
Definition: e_pow.c:98
static const double L2
Definition: e_pow.c:82
static const double zero
Definition: e_pow.c:74
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
static const double ivln2_l
Definition: e_pow.c:101
static const double huge
Definition: e_pow.c:78
static const double L4
Definition: e_pow.c:84
GLdouble n
static const double cp
Definition: e_pow.c:96
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
static const double cp_h
Definition: e_pow.c:97
static const double ivln2
Definition: e_pow.c:99
static const double P5
Definition: e_pow.c:91
libm_hidden_def(scalbln)
Definition: s_scalbn.c:66
static const double P2
Definition: e_pow.c:88
double fabs(double x)
Definition: s_fabs.c:22
GLdouble GLdouble t
Definition: SDL_opengl.h:2071