4.3BSD/usr/contrib/icon/rt/log.c

/*
 * log - calculate the natural logarithm of arg.
 *
 * The coefficients are #2705 from Hart & Cheney. (19.38D)
 */

#include <errno.h>
#include <math.h>

int   errno;
double   frexp();
static double   log2   = 0.693147180559945309e0;
static double   ln10   = 2.302585092994045684;
static double   sqrto2   = 0.707106781186547524e0;
static double   p0   = -.240139179559210510e2;
static double   p1   = 0.309572928215376501e2;
static double   p2   = -.963769093368686593e1;
static double   p3   = 0.421087371217979714e0;
static double   q0   = -.120069589779605255e2;
static double   q1   = 0.194809660700889731e2;
static double   q2   = -.891110902798312337e1;

double
log(arg)
double arg;
   {
   double x, z, zsq, temp;
   int exp;

   if(arg <= 0.) {
      errno = EDOM;
      return(-HUGE);
   }
   x = frexp(arg,&exp);
   while(x<0.5) {
      x = x*2;
      exp = exp-1;
   }
   if(x<sqrto2) {
      x = 2*x;
      exp = exp-1;
   }

   z = (x-1)/(x+1);
   zsq = z*z;

   temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
   temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0);
   temp = temp*z + exp*log2;
   return(temp);
   }

double
log10(arg)
double arg;
   {

   return(log(arg)/ln10);
   }