This is exp.c in view mode; [Download] [Up]
#include <stdio.h>
static double log2 = 0.6931471805599453094;
static double ilog2 = 1./0.6931471805599453094;
double jv_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 ) {
fprintf(stderr,"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.