1*37da2899SCharles.Forsyth /* 2*37da2899SCharles.Forsyth exp returns the exponential function of its 3*37da2899SCharles.Forsyth floating-point argument. 4*37da2899SCharles.Forsyth 5*37da2899SCharles.Forsyth The coefficients are #1069 from Hart and Cheney. (22.35D) 6*37da2899SCharles.Forsyth */ 7*37da2899SCharles.Forsyth 8*37da2899SCharles.Forsyth #include <lib9.h> 9*37da2899SCharles.Forsyth 10*37da2899SCharles.Forsyth #define p0 .2080384346694663001443843411e7 11*37da2899SCharles.Forsyth #define p1 .3028697169744036299076048876e5 12*37da2899SCharles.Forsyth #define p2 .6061485330061080841615584556e2 13*37da2899SCharles.Forsyth #define q0 .6002720360238832528230907598e7 14*37da2899SCharles.Forsyth #define q1 .3277251518082914423057964422e6 15*37da2899SCharles.Forsyth #define q2 .1749287689093076403844945335e4 16*37da2899SCharles.Forsyth #define log2e 1.4426950408889634073599247 17*37da2899SCharles.Forsyth #define sqrt2 1.4142135623730950488016887 18*37da2899SCharles.Forsyth #define maxf 10000 19*37da2899SCharles.Forsyth 20*37da2899SCharles.Forsyth double exp(double arg)21*37da2899SCharles.Forsythexp(double arg) 22*37da2899SCharles.Forsyth { 23*37da2899SCharles.Forsyth double fract, temp1, temp2, xsq; 24*37da2899SCharles.Forsyth int ent; 25*37da2899SCharles.Forsyth 26*37da2899SCharles.Forsyth if(arg == 0) 27*37da2899SCharles.Forsyth return 1; 28*37da2899SCharles.Forsyth if(arg < -maxf) 29*37da2899SCharles.Forsyth return 0; 30*37da2899SCharles.Forsyth if(arg > maxf) 31*37da2899SCharles.Forsyth return Inf(1); 32*37da2899SCharles.Forsyth arg *= log2e; 33*37da2899SCharles.Forsyth ent = floor(arg); 34*37da2899SCharles.Forsyth fract = (arg-ent) - 0.5; 35*37da2899SCharles.Forsyth xsq = fract*fract; 36*37da2899SCharles.Forsyth temp1 = ((p2*xsq+p1)*xsq+p0)*fract; 37*37da2899SCharles.Forsyth temp2 = ((xsq+q2)*xsq+q1)*xsq + q0; 38*37da2899SCharles.Forsyth return ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent); 39*37da2899SCharles.Forsyth } 40