SDL  2.0
e_exp.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_exp(x)
13  * Returns the exponential of x.
14  *
15  * Method
16  * 1. Argument reduction:
17  * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
18  * Given x, find r and integer k such that
19  *
20  * x = k*ln2 + r, |r| <= 0.5*ln2.
21  *
22  * Here r will be represented as r = hi-lo for better
23  * accuracy.
24  *
25  * 2. Approximation of exp(r) by a special rational function on
26  * the interval [0,0.34658]:
27  * Write
28  * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
29  * We use a special Reme algorithm on [0,0.34658] to generate
30  * a polynomial of degree 5 to approximate R. The maximum error
31  * of this polynomial approximation is bounded by 2**-59. In
32  * other words,
33  * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
34  * (where z=r*r, and the values of P1 to P5 are listed below)
35  * and
36  * | 5 | -59
37  * | 2.0+P1*z+...+P5*z - R(z) | <= 2
38  * | |
39  * The computation of exp(r) thus becomes
40  * 2*r
41  * exp(r) = 1 + -------
42  * R - r
43  * r*R1(r)
44  * = 1 + r + ----------- (for better accuracy)
45  * 2 - R1(r)
46  * where
47  * 2 4 10
48  * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
49  *
50  * 3. Scale back to obtain exp(x):
51  * From step 1, we have
52  * exp(x) = 2^k * exp(r)
53  *
54  * Special cases:
55  * exp(INF) is INF, exp(NaN) is NaN;
56  * exp(-INF) is 0, and
57  * for finite argument, only exp(0)=1 is exact.
58  *
59  * Accuracy:
60  * according to an error analysis, the error is always less than
61  * 1 ulp (unit in the last place).
62  *
63  * Misc. info.
64  * For IEEE double
65  * if x > 7.09782712893383973096e+02 then exp(x) overflow
66  * if x < -7.45133219101941108420e+02 then exp(x) underflow
67  *
68  * Constants:
69  * The hexadecimal values are the intended ones for the following
70  * constants. The decimal values may be used, provided that the
71  * compiler will convert from decimal to binary accurately enough
72  * to produce the hexadecimal values shown.
73  */
74 
75 #include "math_libm.h"
76 #include "math_private.h"
77 
78 #ifdef __WATCOMC__ /* Watcom defines huge=__huge */
79 #undef huge
80 #endif
81 
82 static const double
83 one = 1.0,
84 halF[2] = {0.5,-0.5,},
85 huge = 1.0e+300,
86 twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
87 o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
88 u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
89 ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
90  -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
91 ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
92  -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
93 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
94 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
95 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
96 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
97 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
98 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
99 
100 double __ieee754_exp(double x) /* default IEEE double exp */
101 {
102  double y;
103  double hi = 0.0;
104  double lo = 0.0;
105  double c;
106  double t;
107  int32_t k=0;
108  int32_t xsb;
109  u_int32_t hx;
110 
111  GET_HIGH_WORD(hx,x);
112  xsb = (hx>>31)&1; /* sign bit of x */
113  hx &= 0x7fffffff; /* high word of |x| */
114 
115  /* filter out non-finite argument */
116  if(hx >= 0x40862E42) { /* if |x|>=709.78... */
117  if(hx>=0x7ff00000) {
118  u_int32_t lx;
119  GET_LOW_WORD(lx,x);
120  if(((hx&0xfffff)|lx)!=0)
121  return x+x; /* NaN */
122  else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
123  }
124  #if 1
125  if(x > o_threshold) return huge*huge; /* overflow */
126  #else /* !!! FIXME: check this: "huge * huge" is a compiler warning, maybe they wanted +Inf? */
127  if(x > o_threshold) return INFINITY; /* overflow */
128  #endif
129 
130  if(x < u_threshold) return twom1000*twom1000; /* underflow */
131  }
132 
133  /* argument reduction */
134  if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
135  if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
136  hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
137  } else {
138  k = (int32_t) (invln2*x+halF[xsb]);
139  t = k;
140  hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
141  lo = t*ln2LO[0];
142  }
143  x = hi - lo;
144  }
145  else if(hx < 0x3e300000) { /* when |x|<2**-28 */
146  if(huge+x>one) return one+x;/* trigger inexact */
147  }
148  else k = 0;
149 
150  /* x is now in primary range */
151  t = x*x;
152  c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
153  if(k==0) return one-((x*c)/(c-2.0)-x);
154  else y = one-((lo-(x*c)/(2.0-c))-hi);
155  if(k >= -1021) {
156  u_int32_t hy;
157  GET_HIGH_WORD(hy,y);
158  SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */
159  return y;
160  } else {
161  u_int32_t hy;
162  GET_HIGH_WORD(hy,y);
163  SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
164  return y*twom1000;
165  }
166 }
167 
168 /*
169  * wrapper exp(x)
170  */
171 #ifndef _IEEE_LIBM
172 double exp(double x)
173 {
174  static const double o_threshold = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
175  static const double u_threshold = -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */
176 
177  double z = __ieee754_exp(x);
178  if (_LIB_VERSION == _IEEE_)
179  return z;
180  if (isfinite(x)) {
181  if (x > o_threshold)
182  return __kernel_standard(x, x, 6); /* exp overflow */
183  if (x < u_threshold)
184  return __kernel_standard(x, x, 7); /* exp underflow */
185  }
186  return z;
187 }
188 #else
190 #endif
191 libm_hidden_def(exp)
P2
static const double P2
Definition: e_exp.c:95
c
const GLubyte * c
Definition: SDL_opengl_glext.h:11096
P4
static const double P4
Definition: e_exp.c:97
math_private.h
huge
static const double huge
Definition: e_exp.c:85
z
GLdouble GLdouble z
Definition: SDL_opengl_glext.h:407
o_threshold
static const double o_threshold
Definition: e_exp.c:87
libm_hidden_def
libm_hidden_def(scalbln)
Definition: s_scalbn.c:66
SET_HIGH_WORD
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
strong_alias
#define strong_alias(x, y)
Definition: math_private.h:28
math_libm.h
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
u_threshold
static const double u_threshold
Definition: e_exp.c:88
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
t
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
ln2HI
static const double ln2HI[2]
Definition: e_exp.c:89
P1
static const double P1
Definition: e_exp.c:94
y
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
twom1000
static const double twom1000
Definition: e_exp.c:86
P5
static const double P5
Definition: e_exp.c:98
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
P3
static const double P3
Definition: e_exp.c:96
__ieee754_exp
double __ieee754_exp(double x)
Definition: e_exp.c:100
GET_LOW_WORD
#define GET_LOW_WORD(i, d)
Definition: math_private.h:118
invln2
static const double invln2
Definition: e_exp.c:93
halF
static const double halF[2]
Definition: e_exp.c:84
ln2LO
static const double ln2LO[2]
Definition: e_exp.c:91
GET_HIGH_WORD
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
one
static const double one
Definition: e_exp.c:83