1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
/* expf.c
*
* Exponential function
*
*
*
* SYNOPSIS:
*
* float x, y, expf();
*
* y = expf( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A polynomial is used to approximate exp(f)
* in the basic range [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +- MAXLOG 100000 1.7e-7 2.8e-8
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* expf underflow x < MINLOGF 0.0
* expf overflow x > MAXLOGF MAXNUMF
*
*/
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Single precision exponential function.
* test interval: [-0.5, +0.5]
* trials: 80000
* peak relative error: 7.6e-8
* rms relative error: 2.8e-8
*/
#include <math.h>
extern float LOG2EF, MAXLOGF, MINLOGF, MAXNUMF;
static float C1 = 0.693359375;
static float C2 = -2.12194440e-4;
float floorf( float ), ldexpf( float, int );
float expf( float xx )
{
float x, z;
int n;
x = xx;
if( x > MAXLOGF)
{
mtherr( "expf", OVERFLOW );
return( MAXNUMF );
}
if( x < MINLOGF )
{
mtherr( "expf", UNDERFLOW );
return(0.0);
}
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
z = floorf( LOG2EF * x + 0.5 ); /* floor() truncates toward -infinity. */
x -= z * C1;
x -= z * C2;
n = z;
z = x * x;
/* Theoretical peak relative error in [-0.5, +0.5] is 4.2e-9. */
z =
((((( 1.9875691500E-4 * x
+ 1.3981999507E-3) * x
+ 8.3334519073E-3) * x
+ 4.1665795894E-2) * x
+ 1.6666665459E-1) * x
+ 5.0000001201E-1) * z
+ x
+ 1.0;
/* multiply by power of 2 */
x = ldexpf( z, n );
return( x );
}
|
