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