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
82static const double
83one = 1.0,
84halF[2] = {0.5,-0.5,},
85huge = 1.0e+300,
86twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
87o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
88u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
89ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
90 -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
91ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
92 -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
93invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
94P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
95P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
96P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
97P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
98P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
99
100double __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
172double 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
signed int int32_t
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
const GLubyte * c
GLdouble GLdouble z
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 halF[2]
Definition: e_exp.c:84
static const double invln2
Definition: e_exp.c:93
static const double P1
Definition: e_exp.c:94
static const double twom1000
Definition: e_exp.c:86
static const double ln2LO[2]
Definition: e_exp.c:91
static const double ln2HI[2]
Definition: e_exp.c:89
static const double u_threshold
Definition: e_exp.c:88
static const double P3
Definition: e_exp.c:96
static const double o_threshold
Definition: e_exp.c:87
static const double huge
Definition: e_exp.c:85
static const double P2
Definition: e_exp.c:95
static const double one
Definition: e_exp.c:83
double __ieee754_exp(double x)
Definition: e_exp.c:100
static const double P4
Definition: e_exp.c:97
static const double P5
Definition: e_exp.c:98
unsigned int u_int32_t
Definition: math_private.h:31
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
#define strong_alias(x, y)
Definition: math_private.h:28
#define GET_LOW_WORD(i, d)
Definition: math_private.h:118
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
libm_hidden_def(scalbln)
Definition: s_scalbn.c:66