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 <math.h> 9 #include <errno.h> 10 11 #define p0 .2080384346694663001443843411e7 12 #define p1 .3028697169744036299076048876e5 13 #define p2 .6061485330061080841615584556e2 14 #define q0 .6002720360238832528230907598e7 15 #define q1 .3277251518082914423057964422e6 16 #define q2 .1749287689093076403844945335e4 17 #define log2e 1.4426950408889634073599247 18 #define sqrt2 1.4142135623730950488016887 19 #define maxf 10000 20 21 double exp(double arg)22exp(double arg) 23 { 24 double fract, temp1, temp2, xsq; 25 int ent; 26 27 if(arg == 0) 28 return 1; 29 if(arg < -maxf) { 30 errno = ERANGE; 31 return 0; 32 } 33 if(arg > maxf) { 34 errno = ERANGE; 35 return HUGE_VAL; 36 } 37 arg *= log2e; 38 ent = floor(arg); 39 fract = (arg-ent) - 0.5; 40 xsq = fract*fract; 41 temp1 = ((p2*xsq+p1)*xsq+p0)*fract; 42 temp2 = ((xsq+q2)*xsq+q1)*xsq + q0; 43 return ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent); 44 } 45