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 | } |