[444] | 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 */ |
---|