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
123
124
125
126
127
128
129
130
131
132
133
|
/* @(#)z_exp.c 1.0 98/08/13 */
/******************************************************************
* The following routines are coded directly from the algorithms
* and coefficients given in "Software Manual for the Elementary
* Functions" by William J. Cody, Jr. and William Waite, Prentice
* Hall, 1980.
******************************************************************/
/*
FUNCTION
<<exp>>, <<expf>>---exponential
INDEX
exp
INDEX
expf
ANSI_SYNOPSIS
#include <math.h>
double exp(double <[x]>);
float expf(float <[x]>);
TRAD_SYNOPSIS
#include <math.h>
double exp(<[x]>);
double <[x]>;
float expf(<[x]>);
float <[x]>;
DESCRIPTION
<<exp>> and <<expf>> calculate the exponential of <[x]>, that is,
@ifnottex
e raised to the power <[x]> (where e
@end ifnottex
@tex
$e^x$ (where $e$
@end tex
is the base of the natural system of logarithms, approximately 2.71828).
RETURNS
On success, <<exp>> and <<expf>> return the calculated value.
If the result underflows, the returned value is <<0>>. If the
result overflows, the returned value is <<HUGE_VAL>>. In
either case, <<errno>> is set to <<ERANGE>>.
PORTABILITY
<<exp>> is ANSI C. <<expf>> is an extension.
*/
/*****************************************************************
* Exponential Function
*
* Input:
* x - floating point value
*
* Output:
* e raised to x.
*
* Description:
* This routine returns e raised to the xth power.
*
*****************************************************************/
#include <float.h>
#include "fdlibm.h"
#include "zmath.h"
#ifndef _DOUBLE_IS_32BITS
static const double INV_LN2 = 1.4426950408889634074;
static const double LN2 = 0.6931471805599453094172321;
static const double p[] = { 0.25, 0.75753180159422776666e-2,
0.31555192765684646356e-4 };
static const double q[] = { 0.5, 0.56817302698551221787e-1,
0.63121894374398504557e-3,
0.75104028399870046114e-6 };
double
_DEFUN (exp, (double),
double x)
{
int N;
double g, z, R, P, Q;
switch (numtest (x))
{
case NAN:
errno = EDOM;
return (x);
case INF:
errno = ERANGE;
if (ispos (x))
return (z_infinity.d);
else
return (0.0);
case 0:
return (1.0);
}
/* Check for out of bounds. */
if (x > BIGX || x < SMALLX)
{
errno = ERANGE;
return (x);
}
/* Check for a value too small to calculate. */
if (-z_rooteps < x && x < z_rooteps)
{
return (1.0);
}
/* Calculate the exponent. */
if (x < 0.0)
N = (int) (x * INV_LN2 - 0.5);
else
N = (int) (x * INV_LN2 + 0.5);
/* Construct the mantissa. */
g = x - N * LN2;
z = g * g;
P = g * ((p[2] * z + p[1]) * z + p[0]);
Q = ((q[3] * z + q[2]) * z + q[1]) * z + q[0];
R = 0.5 + P / (Q - P);
/* Return the floating point value. */
N++;
return (ldexp (R, N));
}
#endif /* _DOUBLE_IS_32BITS */
|
