This is exp.c in view mode; [Download] [Up]
static double log2 = 0.6931471805599453094; static double ilog2 = 1./0.6931471805599453094; double exp(x) double x; { /* routine for the evaluation of the exp (IEEE) Method: 1: y = x/log(2) -> exp(x) = 2^y; 2: z = n+y, -0.5 <= y < 0.5 3: exp(x) = exp(z*log(2))*2^n For the last series we need only terms till 13! Routine made by J.A.M.Vermaseren 20-apr-1992 */ double y,z,scr; long n; int m,sflag; y = x*ilog2; if ( y < 0 ) { sflag = 1; y = -y; } else sflag = 0; scr = y + 0.5; n = *((long *)(&scr)); m = (n >> 20) - 0x3FF; if ( m > 7 ) { printf("Overflow in exponent of %e\n",x); exit(-1); } n &= 0xFFFFF; n |= 0x100000; n >>= (20-m); z = (y - n)*log2; if ( sflag ) { n = -n; z = -z; } scr = 0.99999999999999999693+z*( 0.99999999999999999850+z*( 0.50000000000000183922+z*( 0.16666666666666701804+z*( 0.04166666666648804594+z*( 0.00833333333330984930+z*( 0.00138888889523274996+z*( 0.00019841269908532560+z*( 0.00002480148546912669+z*( 0.00000275572255877329+z*( 0.00000027632644397330+z*( 0.00000002511466084085))))))))))); n <<= 20; *((long *)(&scr)) += n; return(scr); }
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.