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

 libm_hidden_proto (scalbn)
 
double attribute_hidden __ieee754_pow (double x, double y)
 

Variables

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

Function Documentation

◆ __ieee754_pow()

double attribute_hidden __ieee754_pow ( double  x,
double  y 
)

Definition at line 106 of file e_pow.c.

References __ieee754_sqrt, cp, cp_h, cp_l, dp_h, dp_l, EXTRACT_WORDS, fabs, GET_HIGH_WORD, huge_val, 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.

109  {
110  double z, ax, z_h, z_l, p_h, p_l;
111  double y1, t1, t2, r, s, t, u, v, w;
112  int32_t i, j, k, yisint, n;
113  int32_t hx, hy, ix, iy;
114  u_int32_t lx, ly;
115 
116  EXTRACT_WORDS(hx, lx, x);
117  EXTRACT_WORDS(hy, ly, y);
118  ix = hx & 0x7fffffff;
119  iy = hy & 0x7fffffff;
120 
121  /* y==zero: x**0 = 1 */
122  if ((iy | ly) == 0)
123  return one;
124 
125  /* +-NaN return x+y */
126  if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
127  iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
128  return x + y;
129 
130  /* determine if y is an odd int when x < 0
131  * yisint = 0 ... y is not an integer
132  * yisint = 1 ... y is an odd int
133  * yisint = 2 ... y is an even int
134  */
135  yisint = 0;
136  if (hx < 0) {
137  if (iy >= 0x43400000)
138  yisint = 2; /* even integer y */
139  else if (iy >= 0x3ff00000) {
140  k = (iy >> 20) - 0x3ff; /* exponent */
141  if (k > 20) {
142  j = ly >> (52 - k);
143  if ((j << (52 - k)) == ly)
144  yisint = 2 - (j & 1);
145  } else if (ly == 0) {
146  j = iy >> (20 - k);
147  if ((j << (20 - k)) == iy)
148  yisint = 2 - (j & 1);
149  }
150  }
151  }
152 
153  /* special value of y */
154  if (ly == 0) {
155  if (iy == 0x7ff00000) { /* y is +-inf */
156  if (((ix - 0x3ff00000) | lx) == 0)
157  return y - y; /* inf**+-1 is NaN */
158  else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
159  return (hy >= 0) ? y : zero;
160  else /* (|x|<1)**-,+inf = inf,0 */
161  return (hy < 0) ? -y : zero;
162  }
163  if (iy == 0x3ff00000) { /* y is +-1 */
164  if (hy < 0)
165  return one / x;
166  else
167  return x;
168  }
169  if (hy == 0x40000000)
170  return x * x; /* y is 2 */
171  if (hy == 0x3fe00000) { /* y is 0.5 */
172  if (hx >= 0) /* x >= +0 */
173  return __ieee754_sqrt(x);
174  }
175  }
176 
177  ax = fabs(x);
178  /* special value of x */
179  if (lx == 0) {
180  if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
181  z = ax; /* x is +-0,+-inf,+-1 */
182  if (hy < 0)
183  z = one / z; /* z = (1/|x|) */
184  if (hx < 0) {
185  if (((ix - 0x3ff00000) | yisint) == 0) {
186  z = (z - z) / (z - z); /* (-1)**non-int is NaN */
187  } else if (yisint == 1)
188  z = -z; /* (x<0)**odd = -(|x|**odd) */
189  }
190  return z;
191  }
192  }
193 
194  /* (x<0)**(non-int) is NaN */
195  if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
196  return (x - x) / (x - x);
197 
198  /* |y| is huge */
199  if (iy > 0x41e00000) { /* if |y| > 2**31 */
200  if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */
201  if (ix <= 0x3fefffff)
202  return (hy < 0) ? huge_val * huge_val : tiny * tiny;
203  if (ix >= 0x3ff00000)
204  return (hy > 0) ? huge_val * huge_val : tiny * tiny;
205  }
206  /* over/underflow if x is not close to one */
207  if (ix < 0x3fefffff)
208  return (hy < 0) ? huge_val * huge_val : tiny * tiny;
209  if (ix > 0x3ff00000)
210  return (hy > 0) ? huge_val * huge_val : tiny * tiny;
211  /* now |1-x| is tiny <= 2**-20, suffice to compute
212  log(x) by x-x^2/2+x^3/3-x^4/4 */
213  t = x - 1; /* t has 20 trailing zeros */
214  w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
215  u = ivln2_h * t; /* ivln2_h has 21 sig. bits */
216  v = t * ivln2_l - w * ivln2;
217  t1 = u + v;
218  SET_LOW_WORD(t1, 0);
219  t2 = v - (t1 - u);
220  } else {
221  double s2, s_h, s_l, t_h, t_l;
222  n = 0;
223  /* take care subnormal number */
224  if (ix < 0x00100000) {
225  ax *= two53;
226  n -= 53;
227  GET_HIGH_WORD(ix, ax);
228  }
229  n += ((ix) >> 20) - 0x3ff;
230  j = ix & 0x000fffff;
231  /* determine interval */
232  ix = j | 0x3ff00000; /* normalize ix */
233  if (j <= 0x3988E)
234  k = 0; /* |x|<sqrt(3/2) */
235  else if (j < 0xBB67A)
236  k = 1; /* |x|<sqrt(3) */
237  else {
238  k = 0;
239  n += 1;
240  ix -= 0x00100000;
241  }
242  SET_HIGH_WORD(ax, ix);
243 
244  /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
245  u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
246  v = one / (ax + bp[k]);
247  s = u * v;
248  s_h = s;
249  SET_LOW_WORD(s_h, 0);
250  /* t_h=ax+bp[k] High */
251  t_h = zero;
252  SET_HIGH_WORD(t_h,
253  ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
254  t_l = ax - (t_h - bp[k]);
255  s_l = v * ((u - s_h * t_h) - s_h * t_l);
256  /* compute log(ax) */
257  s2 = s * s;
258  r = s2 * s2 * (L1 +
259  s2 * (L2 +
260  s2 * (L3 +
261  s2 * (L4 + s2 * (L5 + s2 * L6)))));
262  r += s_l * (s_h + s);
263  s2 = s_h * s_h;
264  t_h = 3.0 + s2 + r;
265  SET_LOW_WORD(t_h, 0);
266  t_l = r - ((t_h - 3.0) - s2);
267  /* u+v = s*(1+...) */
268  u = s_h * t_h;
269  v = s_l * t_h + t_l * s;
270  /* 2/(3log2)*(s+...) */
271  p_h = u + v;
272  SET_LOW_WORD(p_h, 0);
273  p_l = v - (p_h - u);
274  z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
275  z_l = cp_l * p_h + p_l * cp + dp_l[k];
276  /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
277  t = (double) n;
278  t1 = (((z_h + z_l) + dp_h[k]) + t);
279  SET_LOW_WORD(t1, 0);
280  t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
281  }
282 
283  s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
284  if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
285  s = -one; /* (-ve)**(odd int) */
286 
287  /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
288  y1 = y;
289  SET_LOW_WORD(y1, 0);
290  p_l = (y - y1) * t1 + y * t2;
291  p_h = y1 * t1;
292  z = p_l + p_h;
293  EXTRACT_WORDS(j, i, z);
294  if (j >= 0x40900000) { /* z >= 1024 */
295  if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
296  return s * huge_val * huge_val; /* overflow */
297  else {
298  if (p_l + ovt > z - p_h)
299  return s * huge_val * huge_val; /* overflow */
300  }
301  } else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */
302  if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */
303  return s * tiny * tiny; /* underflow */
304  else {
305  if (p_l <= z - p_h)
306  return s * tiny * tiny; /* underflow */
307  }
308  }
309  /*
310  * compute 2**(p_h+p_l)
311  */
312  i = j & 0x7fffffff;
313  k = (i >> 20) - 0x3ff;
314  n = 0;
315  if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
316  n = j + (0x00100000 >> (k + 1));
317  k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
318  t = zero;
319  SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
320  n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
321  if (j < 0)
322  n = -n;
323  p_h -= t;
324  }
325  t = p_l + p_h;
326  SET_LOW_WORD(t, 0);
327  u = t * lg2_h;
328  v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
329  z = u + v;
330  w = v - (z - u);
331  t = z * z;
332  t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
333  r = (z * t1) / (t1 - two) - (w + z * w);
334  z = one - (r - z);
335  GET_HIGH_WORD(j, z);
336  j += (n << 20);
337  if ((j >> 20) <= 0)
338  z = scalbn(z, n); /* subnormal output */
339  else
340  SET_HIGH_WORD(z, j);
341  return s * z;
342  }
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:103
ovt
Definition: e_pow.c:95
one
Definition: e_pow.c:78
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
GLdouble GLdouble z
GLdouble s
Definition: SDL_opengl.h:2063
ivln2_l
Definition: e_pow.c:101
const GLdouble * v
Definition: SDL_opengl.h:2064
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
tiny
Definition: e_pow.c:79
two
Definition: e_pow.c:78
cp_l
Definition: e_pow.c:98
P3
Definition: e_pow.c:89
#define __ieee754_sqrt
Definition: math_private.h:42
lg2
Definition: e_pow.c:92
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:131
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
#define scalbn
Definition: math_private.h:40
GLfixed y1
P2
Definition: e_pow.c:88
dp_l[]
Definition: e_pow.c:75
unsigned int u_int32_t
Definition: math_private.h:29
ivln2
Definition: e_pow.c:99
cp
Definition: e_pow.c:96
#define SET_LOW_WORD(d, v)
Definition: math_private.h:141
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:93
L3
Definition: e_pow.c:83
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
GLubyte GLubyte GLubyte GLubyte w
P5
Definition: e_pow.c:91
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
lg2_h
Definition: e_pow.c:93
P4
Definition: e_pow.c:90
#define fabs
Definition: math_private.h:36
L6
Definition: e_pow.c:86
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
ivln2_h
Definition: e_pow.c:100
L5
Definition: e_pow.c:85
lg2_l
Definition: e_pow.c:94
P1
Definition: e_pow.c:87
dp_h[]
Definition: e_pow.c:72
GLdouble n
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
cp_h
Definition: e_pow.c:97
huge_val
Definition: e_pow.c:79
L4
Definition: e_pow.c:84
zero
Definition: e_pow.c:78
L1
Definition: e_pow.c:81
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
L2
Definition: e_pow.c:82
two53
Definition: e_pow.c:78

◆ libm_hidden_proto()

libm_hidden_proto ( scalbn  )

Definition at line 65 of file e_pow.c.

72  { 1.0, 1.5, }, dp_h[] = {
dp_h[]
Definition: e_pow.c:72

Variable Documentation

◆ cp

cp = 9.61796693925975554329e-01

Definition at line 96 of file e_pow.c.

Referenced by __ieee754_pow(), and mmap_resize().

◆ cp_h

cp_h = 9.61796700954437255859e-01

Definition at line 97 of file e_pow.c.

Referenced by __ieee754_pow().

◆ cp_l

cp_l = -7.02846165095275826516e-09

Definition at line 98 of file e_pow.c.

Referenced by __ieee754_pow().

◆ dp_h

dp_h[]
Initial value:
= {
0.0, 5.84962487220764160156e-01,}

Definition at line 72 of file e_pow.c.

Referenced by __ieee754_pow().

◆ dp_l

dp_l[]
Initial value:
= {
0.0, 1.35003920212974897128e-08,}

Definition at line 75 of file e_pow.c.

Referenced by __ieee754_pow().

◆ huge_val

huge_val = 1.0e300

Definition at line 79 of file e_pow.c.

Referenced by __ieee754_pow(), and libm_hidden_proto().

◆ ivln2

ivln2 = 1.44269504088896338700e+00

Definition at line 99 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ivln2_h

ivln2_h = 1.44269502162933349609e+00

Definition at line 100 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ivln2_l

ivln2_l = 1.92596299112661746887e-08

Definition at line 101 of file e_pow.c.

Referenced by __ieee754_pow().

◆ L1

L1 = 5.99999999999994648725e-01

Definition at line 81 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ L2

L2 = 4.28571428578550184252e-01

Definition at line 82 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ L3

L3 = 3.33333329818377432918e-01

Definition at line 83 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ L4

L4 = 2.72728123808534006489e-01

Definition at line 84 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ L5

L5 = 2.30660745775561754067e-01

Definition at line 85 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ L6

L6 = 2.06975017800338417784e-01

Definition at line 86 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

◆ lg2

lg2 = 6.93147180559945286227e-01

Definition at line 92 of file e_pow.c.

Referenced by __ieee754_pow().

◆ lg2_h

lg2_h = 6.93147182464599609375e-01

Definition at line 93 of file e_pow.c.

Referenced by __ieee754_pow().

◆ lg2_l

lg2_l = -1.90465429995776804525e-09

Definition at line 94 of file e_pow.c.

Referenced by __ieee754_pow().

◆ one

one = 1.0

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow().

◆ ovt

ovt = 8.0085662595372944372e-0017

Definition at line 95 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P1

P1 = 1.66666666666666019037e-01

Definition at line 87 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P2

P2 = -2.77777777770155933842e-03

Definition at line 88 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P3

P3 = 6.61375632143793436117e-05

Definition at line 89 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P4

P4 = -1.65339022054652515390e-06

Definition at line 90 of file e_pow.c.

Referenced by __ieee754_pow().

◆ P5

P5 = 4.13813679705723846039e-08

Definition at line 91 of file e_pow.c.

Referenced by __ieee754_pow().

◆ tiny

tiny = 1.0e-300

Definition at line 79 of file e_pow.c.

Referenced by __ieee754_pow(), and libm_hidden_proto().

◆ two

two = 2.0

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow().

◆ two53

two53 = 9007199254740992.0

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow().

◆ zero

zero = 0.0

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_Convert_F32_to_S32_Scalar().