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