SDL  2.0
e_pow.c File Reference
#include "math_libm.h"
#include "math_private.h"
+ Include dependency graph for e_pow.c:

Go to the source code of this file.

Functions

double attribute_hidden __ieee754_pow (double x, double y)
 

Variables

static const double bp [] = {1.0, 1.5,}
 
static const double dp_h [] = { 0.0, 5.84962487220764160156e-01,}
 
static const double dp_l [] = { 0.0, 1.35003920212974897128e-08,}
 
static const double zero = 0.0
 
static const double one = 1.0
 
static const double two = 2.0
 
static const double two53 = 9007199254740992.0
 
static const double huge = 1.0e300
 
static const double tiny = 1.0e-300
 
static const double L1 = 5.99999999999994648725e-01
 
static const double L2 = 4.28571428578550184252e-01
 
static const double L3 = 3.33333329818377432918e-01
 
static const double L4 = 2.72728123808534006489e-01
 
static const double L5 = 2.30660745775561754067e-01
 
static const double L6 = 2.06975017800338417784e-01
 
static const double P1 = 1.66666666666666019037e-01
 
static const double P2 = -2.77777777770155933842e-03
 
static const double P3 = 6.61375632143793436117e-05
 
static const double P4 = -1.65339022054652515390e-06
 
static const double P5 = 4.13813679705723846039e-08
 
static const double lg2 = 6.93147180559945286227e-01
 
static const double lg2_h = 6.93147182464599609375e-01
 
static const double lg2_l = -1.90465429995776804525e-09
 
static const double ovt = 8.0085662595372944372e-0017
 
static const double cp = 9.61796693925975554329e-01
 
static const double cp_h = 9.61796700954437255859e-01
 
static const double cp_l = -7.02846165095275826516e-09
 
static const double ivln2 = 1.44269504088896338700e+00
 
static const double ivln2_h = 1.44269502162933349609e+00
 
static const double ivln2_l = 1.92596299112661746887e-08
 

Function Documentation

◆ __ieee754_pow()

double attribute_hidden __ieee754_pow ( double  x,
double  y 
)

Definition at line 103 of file e_pow.c.

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 }

References __ieee754_sqrt, bp, cp, cp_h, cp_l, dp_h, dp_l, EXTRACT_WORDS, fabs(), GET_HIGH_WORD, huge, i, ivln2, ivln2_h, ivln2_l, j, k, L1, L2, L3, L4, L5, L6, lg2, lg2_h, lg2_l, one, ovt, P1, P2, P3, P4, P5, scalbn, SET_HIGH_WORD, SET_LOW_WORD, tiny, two, two53, and zero.

Variable Documentation

◆ bp

const double bp[] = {1.0, 1.5,}
static

Definition at line 71 of file e_pow.c.

Referenced by __ieee754_pow().

◆ cp

const double cp = 9.61796693925975554329e-01
static

Definition at line 96 of file e_pow.c.

Referenced by __ieee754_pow(), mmap_resize(), and validate_hex().

◆ cp_h

const double cp_h = 9.61796700954437255859e-01
static

Definition at line 97 of file e_pow.c.

Referenced by __ieee754_pow().

◆ cp_l

const double cp_l = -7.02846165095275826516e-09
static

Definition at line 98 of file e_pow.c.

Referenced by __ieee754_pow().

◆ dp_h

const double dp_h[] = { 0.0, 5.84962487220764160156e-01,}
static

Definition at line 72 of file e_pow.c.

Referenced by __ieee754_pow().

◆ dp_l

const double dp_l[] = { 0.0, 1.35003920212974897128e-08,}
static

Definition at line 73 of file e_pow.c.

Referenced by __ieee754_pow().

◆ huge

const double huge = 1.0e300
static

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ivln2

const double ivln2 = 1.44269504088896338700e+00
static

Definition at line 99 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ivln2_h

const double ivln2_h = 1.44269502162933349609e+00
static

Definition at line 100 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ivln2_l

const double ivln2_l = 1.92596299112661746887e-08
static

Definition at line 101 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L1

const double L1 = 5.99999999999994648725e-01
static

Definition at line 81 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L2

const double L2 = 4.28571428578550184252e-01
static

Definition at line 82 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L3

const double L3 = 3.33333329818377432918e-01
static

Definition at line 83 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L4

const double L4 = 2.72728123808534006489e-01
static

Definition at line 84 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L5

const double L5 = 2.30660745775561754067e-01
static

Definition at line 85 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L6

const double L6 = 2.06975017800338417784e-01
static

Definition at line 86 of file e_pow.c.

Referenced by __ieee754_pow().

◆ lg2

const double lg2 = 6.93147180559945286227e-01
static

Definition at line 92 of file e_pow.c.

Referenced by __ieee754_pow().

◆ lg2_h

const double lg2_h = 6.93147182464599609375e-01
static

Definition at line 93 of file e_pow.c.

Referenced by __ieee754_pow().

◆ lg2_l

const double lg2_l = -1.90465429995776804525e-09
static

Definition at line 94 of file e_pow.c.

Referenced by __ieee754_pow().

◆ one

const double one = 1.0
static

Definition at line 75 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ovt

const double ovt = 8.0085662595372944372e-0017
static

Definition at line 95 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P1

const double P1 = 1.66666666666666019037e-01
static

Definition at line 87 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P2

const double P2 = -2.77777777770155933842e-03
static

Definition at line 88 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P3

const double P3 = 6.61375632143793436117e-05
static

Definition at line 89 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P4

const double P4 = -1.65339022054652515390e-06
static

Definition at line 90 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P5

const double P5 = 4.13813679705723846039e-08
static

Definition at line 91 of file e_pow.c.

Referenced by __ieee754_pow().

◆ tiny

const double tiny = 1.0e-300
static

Definition at line 79 of file e_pow.c.

Referenced by __ieee754_pow().

◆ two

const double two = 2.0
static

Definition at line 76 of file e_pow.c.

Referenced by __ieee754_pow().

◆ two53

const double two53 = 9007199254740992.0
static

Definition at line 77 of file e_pow.c.

Referenced by __ieee754_pow().

◆ zero

const double zero = 0.0
static

Definition at line 74 of file e_pow.c.

Referenced by __ieee754_pow().

bp
static const double bp[]
Definition: e_pow.c:71
ivln2_l
static const double ivln2_l
Definition: e_pow.c:101
t1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: SDL_opengl_glext.h:8586
fabs
double fabs(double x)
Definition: s_fabs.c:22
two
static const double two
Definition: e_pow.c:76
cp
static const double cp
Definition: e_pow.c:96
L5
static const double L5
Definition: e_pow.c:85
r
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
P5
static const double P5
Definition: e_pow.c:91
z
GLdouble GLdouble z
Definition: SDL_opengl_glext.h:407
v
const GLdouble * v
Definition: SDL_opengl.h:2064
__ieee754_sqrt
#define __ieee754_sqrt
Definition: math_private.h:48
lg2_h
static const double lg2_h
Definition: e_pow.c:93
zero
static const double zero
Definition: e_pow.c:74
EXTRACT_WORDS
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:99
n
GLdouble n
Definition: SDL_opengl_glext.h:1955
L4
static const double L4
Definition: e_pow.c:84
L2
static const double L2
Definition: e_pow.c:82
SET_HIGH_WORD
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
scalbn
#define scalbn
Definition: math_private.h:46
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
huge
static const double huge
Definition: e_pow.c:78
SET_LOW_WORD
#define SET_LOW_WORD(d, v)
Definition: math_private.h:147
int32_t
signed int int32_t
Definition: SDL_config_windows.h:62
u_int32_t
unsigned int u_int32_t
Definition: math_private.h:31
cp_h
static const double cp_h
Definition: e_pow.c:97
t
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
lg2_l
static const double lg2_l
Definition: e_pow.c:94
ivln2
static const double ivln2
Definition: e_pow.c:99
dp_h
static const double dp_h[]
Definition: e_pow.c:72
P2
static const double P2
Definition: e_pow.c:88
tiny
static const double tiny
Definition: e_pow.c:79
y
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
P4
static const double P4
Definition: e_pow.c:90
k
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 int int int return Display long XVisualInfo int return Display Window Atom long long Bool Atom Atom int unsigned long unsigned long k)
Definition: SDL_x11sym.h:80
dp_l
static const double dp_l[]
Definition: e_pow.c:73
L3
static const double L3
Definition: e_pow.c:83
P1
static const double P1
Definition: e_pow.c:87
s
GLdouble s
Definition: SDL_opengl.h:2063
y1
GLfixed y1
Definition: SDL_opengl_glext.h:4586
two53
static const double two53
Definition: e_pow.c:77
ivln2_h
static const double ivln2_h
Definition: e_pow.c:100
L1
static const double L1
Definition: e_pow.c:81
lg2
static const double lg2
Definition: e_pow.c:92
L6
static const double L6
Definition: e_pow.c:86
GET_HIGH_WORD
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
j
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
one
static const double one
Definition: e_pow.c:75
P3
static const double P3
Definition: e_pow.c:89
cp_l
static const double cp_l
Definition: e_pow.c:98
i
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
ovt
static const double ovt
Definition: e_pow.c:95
w
GLubyte GLubyte GLubyte GLubyte w
Definition: SDL_opengl_glext.h:734