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