quick logarithm

my_log.c

#include <math.h>

double my_log(double x) {
  static const double log2 = 0.693147;

  double y, logy;
  int    n;

  // factor x = y * 2^n, where 1 <= y < 2
  // frexp gives us the range 0.5 <= y < 1
  // so adjust by a factor of 2
  y = frexp(x, &n) * 2;

  // \int_1^y \frac{1}{x} dx = log(y) - log(1) = log(y)
  // by the fundamental theorem of calculus. use simpson's
  // rule for a cheap approximation
  #define inv(x) (1.0 / (x))
  #define simp(f, a, b) ((b - a) / 6.0 * (f(a) + 4 * f((a + b) / 2.0) + f(b)))
  logy = simp(inv, 1.0, y);

  // log(x) = log(y * 2^n) = log(y) + n * log(2)
  // n has to be adjusted as well
  return logy + (n - 1) * log2;
}