This is log.c in view mode; [Download] [Up]
#include <stdio.h>
double jv_log(x)
double x;
{
/*
Routine for log(x) in IEEE double precision.
Note that we have to hack around in the exponent!
Steps: 1 x = y*2^n, with y between sqrt(2) and 1/sqrt(2)
2 convert to z = (y-1)/(y+1)
3 ln = 2*z*sum(i=0,9, z^(2*i)/(2*i+1)) + n*ln(2)
Routine made by J.A.M. Vermaseren 20-apr-1992.
*/
double scr, y, z, z2, sum, jv_sqrt();
short n;
if ( x < 0 ) {
fprintf(stderr,"log of a negative number: %d\n",x);
exit(-1);
}
scr = x;
n = ( (*((short *)(&scr))) & 0x7FF0 ) - 0x3FF0;
*((short *)(&scr)) -= n;
n >>= 4;
if ( scr > 1.4142135623730950488 ) {
*((short *)(&scr)) -= 0x10;
n++;
}
y = scr;
z = (y-1.)/(y+1.);
z2 = z*z;
sum = 2.*z*(0.99999999999999999886+z2*(
0.33333333333333825805+z2*(
0.19999999999649522967+z2*(
0.14285714380662481221+z2*(
0.11111098494595009814+z2*(
0.09091817553294820737+z2*(
0.07656220982777783836+z2*(
0.07405280893003482131))))))))
+n*0.6931471805599453094;
return(sum);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.