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