/*
 * Brian Gaeke <brg@eecs.berkeley.edu>
 *
 *
 * Approximation 2503 from Hart & Cheney. pp 182 & 205.
 * 1. ln x = ln 2 lg x
 * 2. x = f 2^n so lg x = n + lg f   -- f = frexp(x, &n); 
 *    (ln x = (ln 2) (n + lg f))
 * 3. lg f, where f in [1/2,1) = P(x)/Q(x) is
 *    -3.51287525 +8.14029839 f +-8.4829424 f^2 +5.16135624 f^3 +-1.30592456 f^4
 *    which is
 *    f * (8.14029839  +-8.4829424 f^1 +5.16135624 f^2 +-1.30592456 f^3) - 3.51287525
 *    = f * (8.14029839  + f*(-8.4829424 +5.16135624 f^1 +-1.30592456 f^2)) - 3.51287525
 *    = f * (8.14029839  + f*(-8.4829424 +f*(5.16135624 +-1.30592456 f^1))) - 3.51287525
 *    = f * (8.14029839  + f*(-8.4829424 +f*(5.16135624 +f*(-1.30592456)))) - 3.51287525
 */
#include <math.h>
double mylog(double x)
{
  const double c0 = -1.30592456;
  const double c1 = 5.16135624;
  const double c2 = 8.4829424;
  const double c3 = 8.14029839;
  const double c4 = 3.51287525;
  const double ln2 = 0.69314718055994530941723212145817657;
  double f, t0, t1, t2, t3, t4, t5;
  int n;

  f = frexp(x, &n);
  t0 = c0*f + c1;
  t1 = t0*f - c2;
  t2 = t1*f + c3;
  t3 = t2*f - c4;
  t4 = n + t3;
  t5 = ln2 * t4;
  return t5;
}