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