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 <u.h> 9 #include <libc.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 return 0; 31 if(arg > maxf) 32 return Inf(1); 33 arg *= log2e; 34 ent = floor(arg); 35 fract = (arg-ent) - 0.5; 36 xsq = fract*fract; 37 temp1 = ((p2*xsq+p1)*xsq+p0)*fract; 38 temp2 = ((xsq+q2)*xsq+q1)*xsq + q0; 39 return ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent); 40 } 41