| 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 */ |
|---|
