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 99 of file e_pow.c.

References __ieee754_pow(), __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, libm_hidden_def(), one, ovt, P1, P2, P3, P4, P5, scalbn, SET_HIGH_WORD, SET_LOW_WORD, strong_alias, tiny, two, two53, and zero.

Referenced by __ieee754_pow().

{
double z,ax,z_h,z_l,p_h,p_l;
double y1,t1,t2,r,s,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
EXTRACT_WORDS(hx,lx,x);
/* x==1: 1**y = 1 (even if y is NaN) */
if (hx==0x3ff00000 && lx==0) {
return x;
}
ix = hx&0x7fffffff;
EXTRACT_WORDS(hy,ly,y);
iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
yisint = 0;
if(hx<0) {
if(iy>=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
if((j<<(52-k))==ly) yisint = 2-(j&1);
} else if(ly==0) {
j = iy>>(20-k);
if((j<<(20-k))==iy) yisint = 2-(j&1);
}
}
}
/* special value of y */
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if (((ix-0x3ff00000)|lx)==0)
return one; /* +-1**+-inf is 1 (yes, weird rule) */
if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return (hy>=0) ? y : zero;
/* (|x|<1)**-,+inf = inf,0 */
return (hy<0) ? -y : zero;
}
if(iy==0x3ff00000) { /* y is +-1 */
if(hy<0) return one/x; else return x;
}
if(hy==0x40000000) return x*x; /* y is 2 */
if(hy==0x3fe00000) { /* y is 0.5 */
if(hx>=0) /* x >= +0 */
return __ieee754_sqrt(x);
}
}
ax = fabs(x);
/* special value of x */
if(lx==0) {
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /* z = (1/|x|) */
if(hx<0) {
if(((ix-0x3ff00000)|yisint)==0) {
z = (z-z)/(z-z); /* (-1)**non-int is NaN */
} else if(yisint==1)
z = -z; /* (x<0)**odd = -(|x|**odd) */
}
return z;
}
}
/* (x<0)**(non-int) is NaN */
if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
/* |y| is huge */
if(iy>0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
/* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */
t = x-1; /* t has 20 trailing zeros */
w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
t1 = u+v;
SET_LOW_WORD(t1,0);
t2 = v-(t1-u);
} else {
double s2,s_h,s_l,t_h,t_l;
n = 0;
/* take care subnormal number */
if(ix<0x00100000)
{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
else {k=0;n+=1;ix -= 0x00100000;}
SET_HIGH_WORD(ax,ix);
/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
v = one/(ax+bp[k]);
s = u*v;
s_h = s;
SET_LOW_WORD(s_h,0);
/* t_h=ax+bp[k] High */
t_h = zero;
SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
s2 = s*s;
r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
r += s_l*(s_h+s);
s2 = s_h*s_h;
t_h = 3.0+s2+r;
SET_LOW_WORD(t_h,0);
t_l = r-((t_h-3.0)-s2);
/* u+v = s*(1+...) */
u = s_h*t_h;
v = s_l*t_h+t_l*s;
/* 2/(3log2)*(s+...) */
p_h = u+v;
SET_LOW_WORD(p_h,0);
p_l = v-(p_h-u);
z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
z_l = cp_l*p_h+p_l*cp+dp_l[k];
/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
t = (double)n;
t1 = (((z_h+z_l)+dp_h[k])+t);
SET_LOW_WORD(t1,0);
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
s = -one;/* (-ve)**(odd int) */
/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
y1 = y;
SET_LOW_WORD(y1,0);
p_l = (y-y1)*t1+y*t2;
p_h = y1*t1;
z = p_l+p_h;
EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */
else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
}
}
/*
* compute 2**(p_h+p_l)
*/
i = j&0x7fffffff;
k = (i>>20)-0x3ff;
n = 0;
if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j+(0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
t = zero;
SET_HIGH_WORD(t,n&~(0x000fffff>>k));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
SET_LOW_WORD(t,0);
u = t*lg2_h;
v = (p_l-(t-p_h))*lg2+t*lg2_l;
z = u+v;
w = v-(z-u);
t = z*z;
t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
r = (z*t1)/(t1-two)-(w+z*w);
z = one-(r-z);
GET_HIGH_WORD(j,z);
j += (n<<20);
if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
else SET_HIGH_WORD(z,j);
return s*z;
}
GET_HIGH_WORD
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:108
r
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
two53
static const double two53
Definition: e_pow.c:73
z
GLdouble GLdouble z
Definition: SDL_opengl_glext.h:404
dp_l
static const double dp_l[]
Definition: e_pow.c:69
s
GLdouble s
Definition: SDL_opengl.h:2063
v
const GLdouble * v
Definition: SDL_opengl.h:2064
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
int32_t
signed int int32_t
Definition: SDL_config_windows.h:62
ivln2_h
static const double ivln2_h
Definition: e_pow.c:96
P3
static const double P3
Definition: e_pow.c:85
L1
static const double L1
Definition: e_pow.c:77
L6
static const double L6
Definition: e_pow.c:82
lg2
static const double lg2
Definition: e_pow.c:88
lg2_l
static const double lg2_l
Definition: e_pow.c:90
ovt
static const double ovt
Definition: e_pow.c:91
tiny
static const double tiny
Definition: e_pow.c:75
bp
static const double bp[]
Definition: e_pow.c:67
__ieee754_sqrt
#define __ieee754_sqrt
Definition: math_private.h:47
SET_HIGH_WORD
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:136
t1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: SDL_opengl_glext.h:8583
scalbn
#define scalbn
Definition: math_private.h:45
dp_h
static const double dp_h[]
Definition: e_pow.c:68
y1
GLfixed y1
Definition: SDL_opengl_glext.h:4583
L3
static const double L3
Definition: e_pow.c:79
P1
static const double P1
Definition: e_pow.c:83
P4
static const double P4
Definition: e_pow.c:86
u_int32_t
unsigned int u_int32_t
Definition: math_private.h:31
L5
static const double L5
Definition: e_pow.c:81
two
static const double two
Definition: e_pow.c:72
SET_LOW_WORD
#define SET_LOW_WORD(d, v)
Definition: math_private.h:146
EXTRACT_WORDS
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:98
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:71
w
GLubyte GLubyte GLubyte GLubyte w
Definition: SDL_opengl_glext.h:731
lg2_h
static const double lg2_h
Definition: e_pow.c:89
y
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
cp_l
static const double cp_l
Definition: e_pow.c:94
L2
static const double L2
Definition: e_pow.c:78
zero
static const double zero
Definition: e_pow.c:70
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
ivln2_l
static const double ivln2_l
Definition: e_pow.c:97
huge
static const double huge
Definition: e_pow.c:74
L4
static const double L4
Definition: e_pow.c:80
n
GLdouble n
Definition: SDL_opengl_glext.h:1952
cp
static const double cp
Definition: e_pow.c:92
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 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
static const double cp_h
Definition: e_pow.c:93
ivln2
static const double ivln2
Definition: e_pow.c:95
P5
static const double P5
Definition: e_pow.c:87
P2
static const double P2
Definition: e_pow.c:84
fabs
double fabs(double x)
Definition: s_fabs.c:22
t
GLdouble GLdouble t
Definition: SDL_opengl.h:2071

## Variable Documentation

## bp

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

static

Definition at line 67 of file e_pow.c.

Referenced by __ieee754_pow().

## cp

const double cp = 9.61796693925975554329e-01

static

Definition at line 92 of file e_pow.c.

Referenced by __ieee754_pow(), and mmap_resize().

## cp_h

const double cp_h = 9.61796700954437255859e-01

static

Definition at line 93 of file e_pow.c.

Referenced by __ieee754_pow().

## cp_l

const double cp_l = -7.02846165095275826516e-09

static

Definition at line 94 of file e_pow.c.

Referenced by __ieee754_pow().

## dp_h

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

static

Definition at line 68 of file e_pow.c.

Referenced by __ieee754_pow().

## dp_l

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

static

Definition at line 69 of file e_pow.c.

Referenced by __ieee754_pow().

## huge

const double huge = 1.0e300

static

Definition at line 74 of file e_pow.c.

Referenced by __ieee754_pow().

## ivln2

const double ivln2 = 1.44269504088896338700e+00

static

Definition at line 95 of file e_pow.c.

Referenced by __ieee754_pow().

## ivln2_h

const double ivln2_h = 1.44269502162933349609e+00

static

Definition at line 96 of file e_pow.c.

Referenced by __ieee754_pow().

## ivln2_l

const double ivln2_l = 1.92596299112661746887e-08

static

Definition at line 97 of file e_pow.c.

Referenced by __ieee754_pow().

## L1

const double L1 = 5.99999999999994648725e-01

static

Definition at line 77 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## L2

const double L2 = 4.28571428578550184252e-01

static

Definition at line 78 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## L3

const double L3 = 3.33333329818377432918e-01

static

Definition at line 79 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## L4

const double L4 = 2.72728123808534006489e-01

static

Definition at line 80 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## L5

const double L5 = 2.30660745775561754067e-01

static

Definition at line 81 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## L6

const double L6 = 2.06975017800338417784e-01

static

Definition at line 82 of file e_pow.c.

Referenced by __ieee754_pow(), and SDL_tolower().

## lg2

const double lg2 = 6.93147180559945286227e-01

static

Definition at line 88 of file e_pow.c.

Referenced by __ieee754_pow().

## lg2_h

const double lg2_h = 6.93147182464599609375e-01

static

Definition at line 89 of file e_pow.c.

Referenced by __ieee754_pow().

## lg2_l

const double lg2_l = -1.90465429995776804525e-09

static

Definition at line 90 of file e_pow.c.

Referenced by __ieee754_pow().

## one

const double one = 1.0

static

Definition at line 71 of file e_pow.c.

Referenced by __ieee754_pow().

## ovt

const double ovt = 8.0085662595372944372e-0017

static

Definition at line 91 of file e_pow.c.

Referenced by __ieee754_pow().

## P1

const double P1 = 1.66666666666666019037e-01

static

Definition at line 83 of file e_pow.c.

Referenced by __ieee754_pow().

## P2

const double P2 = -2.77777777770155933842e-03

static

Definition at line 84 of file e_pow.c.

Referenced by __ieee754_pow().

## P3

const double P3 = 6.61375632143793436117e-05

static

Definition at line 85 of file e_pow.c.

Referenced by __ieee754_pow().

## P4

const double P4 = -1.65339022054652515390e-06

static

Definition at line 86 of file e_pow.c.

Referenced by __ieee754_pow().

## P5

const double P5 = 4.13813679705723846039e-08

static

Definition at line 87 of file e_pow.c.

Referenced by __ieee754_pow().

## tiny

const double tiny = 1.0e-300

static

Definition at line 75 of file e_pow.c.

Referenced by __ieee754_pow().

## two

const double two = 2.0

static

Definition at line 72 of file e_pow.c.

Referenced by __ieee754_pow().

## two53

const double two53 = 9007199254740992.0

static

Definition at line 73 of file e_pow.c.

Referenced by __ieee754_pow().

## zero

const double zero = 0.0

static

Definition at line 70 of file e_pow.c.

Referenced by __ieee754_pow().