projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 94759b7b TL |
1 | /* |
| 2 | exp returns the exponential function of its | |
| 3 | floating-point argument. | |
| 4 | ||
| 5 | The coefficients are #1069 from Hart and Cheney. (22.35D) | |
| 6 | */ | |
| 7 | ||
| 8 | #include <errno.h> | |
| 9 | #include <math.h> | |
| 10 | ||
| 11 | int errno; | |
| 12 | static double p0 = .2080384346694663001443843411e7; | |
| 13 | static double p1 = .3028697169744036299076048876e5; | |
| 14 | static double p2 = .6061485330061080841615584556e2; | |
| 15 | static double q0 = .6002720360238832528230907598e7; | |
| 16 | static double q1 = .3277251518082914423057964422e6; | |
| 17 | static double q2 = .1749287689093076403844945335e4; | |
| 18 | static double log2e = 1.4426950408889634073599247; | |
| 19 | static double sqrt2 = 1.4142135623730950488016887; | |
| 20 | static double maxf = 10000; | |
| 21 | ||
| 22 | double | |
| 23 | exp(arg) | |
| 24 | double arg; | |
| 25 | { | |
| 26 | double fract; | |
| 27 | double temp1, temp2, xsq; | |
| 28 | int ent; | |
| 29 | ||
| 30 | if(arg == 0.) | |
| 31 | return(1.); | |
| 32 | if(arg < -maxf) | |
| 33 | return(0.); | |
| 34 | if(arg > maxf) { | |
| 35 | errno = ERANGE; | |
| 36 | return(HUGE); | |
| 37 | } | |
| 38 | arg *= log2e; | |
| 39 | ent = floor(arg); | |
| 40 | fract = (arg-ent) - 0.5; | |
| 41 | xsq = fract*fract; | |
| 42 | temp1 = ((p2*xsq+p1)*xsq+p0)*fract; | |
| 43 | temp2 = ((1.0*xsq+q2)*xsq+q1)*xsq + q0; | |
| 44 | return(ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent)); | |
| 45 | } |