|
|
1.1 ! root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.