1 | |
---|
2 | /* @(#)z_exp.c 1.0 98/08/13 */ |
---|
3 | /****************************************************************** |
---|
4 | * The following routines are coded directly from the algorithms |
---|
5 | * and coefficients given in "Software Manual for the Elementary |
---|
6 | * Functions" by William J. Cody, Jr. and William Waite, Prentice |
---|
7 | * Hall, 1980. |
---|
8 | ******************************************************************/ |
---|
9 | |
---|
10 | /* |
---|
11 | FUNCTION |
---|
12 | <<exp>>, <<expf>>---exponential |
---|
13 | INDEX |
---|
14 | exp |
---|
15 | INDEX |
---|
16 | expf |
---|
17 | |
---|
18 | SYNOPSIS |
---|
19 | #include <math.h> |
---|
20 | double exp(double <[x]>); |
---|
21 | float expf(float <[x]>); |
---|
22 | |
---|
23 | DESCRIPTION |
---|
24 | <<exp>> and <<expf>> calculate the exponential of <[x]>, that is, |
---|
25 | @ifnottex |
---|
26 | e raised to the power <[x]> (where e |
---|
27 | @end ifnottex |
---|
28 | @tex |
---|
29 | $e^x$ (where $e$ |
---|
30 | @end tex |
---|
31 | is the base of the natural system of logarithms, approximately 2.71828). |
---|
32 | |
---|
33 | RETURNS |
---|
34 | On success, <<exp>> and <<expf>> return the calculated value. |
---|
35 | If the result underflows, the returned value is <<0>>. If the |
---|
36 | result overflows, the returned value is <<HUGE_VAL>>. In |
---|
37 | either case, <<errno>> is set to <<ERANGE>>. |
---|
38 | |
---|
39 | PORTABILITY |
---|
40 | <<exp>> is ANSI C. <<expf>> is an extension. |
---|
41 | |
---|
42 | */ |
---|
43 | |
---|
44 | /***************************************************************** |
---|
45 | * Exponential Function |
---|
46 | * |
---|
47 | * Input: |
---|
48 | * x - floating point value |
---|
49 | * |
---|
50 | * Output: |
---|
51 | * e raised to x. |
---|
52 | * |
---|
53 | * Description: |
---|
54 | * This routine returns e raised to the xth power. |
---|
55 | * |
---|
56 | *****************************************************************/ |
---|
57 | |
---|
58 | #include <float.h> |
---|
59 | #include "fdlibm.h" |
---|
60 | #include "zmath.h" |
---|
61 | |
---|
62 | #ifndef _DOUBLE_IS_32BITS |
---|
63 | |
---|
64 | static const double INV_LN2 = 1.4426950408889634074; |
---|
65 | static const double LN2 = 0.6931471805599453094172321; |
---|
66 | static const double p[] = { 0.25, 0.75753180159422776666e-2, |
---|
67 | 0.31555192765684646356e-4 }; |
---|
68 | static const double q[] = { 0.5, 0.56817302698551221787e-1, |
---|
69 | 0.63121894374398504557e-3, |
---|
70 | 0.75104028399870046114e-6 }; |
---|
71 | |
---|
72 | double |
---|
73 | exp (double x) |
---|
74 | { |
---|
75 | int N; |
---|
76 | double g, z, R, P, Q; |
---|
77 | |
---|
78 | switch (numtest (x)) |
---|
79 | { |
---|
80 | case NAN: |
---|
81 | errno = EDOM; |
---|
82 | return (x); |
---|
83 | case INF: |
---|
84 | errno = ERANGE; |
---|
85 | if (ispos (x)) |
---|
86 | return (z_infinity.d); |
---|
87 | else |
---|
88 | return (0.0); |
---|
89 | case 0: |
---|
90 | return (1.0); |
---|
91 | } |
---|
92 | |
---|
93 | /* Check for out of bounds. */ |
---|
94 | if (x > BIGX || x < SMALLX) |
---|
95 | { |
---|
96 | errno = ERANGE; |
---|
97 | return (x); |
---|
98 | } |
---|
99 | |
---|
100 | /* Check for a value too small to calculate. */ |
---|
101 | if (-z_rooteps < x && x < z_rooteps) |
---|
102 | { |
---|
103 | return (1.0); |
---|
104 | } |
---|
105 | |
---|
106 | /* Calculate the exponent. */ |
---|
107 | if (x < 0.0) |
---|
108 | N = (int) (x * INV_LN2 - 0.5); |
---|
109 | else |
---|
110 | N = (int) (x * INV_LN2 + 0.5); |
---|
111 | |
---|
112 | /* Construct the mantissa. */ |
---|
113 | g = x - N * LN2; |
---|
114 | z = g * g; |
---|
115 | P = g * ((p[2] * z + p[1]) * z + p[0]); |
---|
116 | Q = ((q[3] * z + q[2]) * z + q[1]) * z + q[0]; |
---|
117 | R = 0.5 + P / (Q - P); |
---|
118 | |
---|
119 | /* Return the floating point value. */ |
---|
120 | N++; |
---|
121 | return (ldexp (R, N)); |
---|
122 | } |
---|
123 | |
---|
124 | #endif /* _DOUBLE_IS_32BITS */ |
---|