/*
 * See http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
 *
 * All of these rely on being on a little endian machine, such as an Intel box.
 *
 * These can be _quite_ inaccurate. ~20% in many cases, but being much faster (~7x) may
 * permit more loop iterations of tuning algorithms that only need approximate powers.
 *
 * This version of Ankerl's algorithm has been extended to provide optionally conservative (lower) bounds
 * and also to generate a full linear interpolation across the entire significand rather than 'stair-step'
 * at the expense of performing a 64 bit operation rather than a 32 bit one. This is cheap these days.
 *
 * 'exp' is further improved by using a suggestion by Nic Schraudolph:
 *
 * "You can get a much better approximation (piecewise rational instead of linear) at
 * the cost of a single floating-point division by using better_exp(x) = exp(x/2)/exp(-x/2),
 * where exp() is my published approximation but you don't need the additive constant anymore,
 * you can use c=0. On machines with hardware division this is very attractive." -- Nic Schraudolph
 *
 * --Edward Kmett
 *
 * TODO: Incorporate the techniques from https://code.google.com/p/fastapprox/ to enable us
 * to calculate more interesting approximate functions. They might need to be generalized to work on
 * Double values where appropriate I suppose.
 *
 * Magic numbers:
 * float /int      : round(1<<23/log(2)) = 12102203,          127<<23 = 1065353216
 * double/int      : round(1<<20/log(2)) = 1512775,          1023<<20 = 1072693248
 * double/long long: round(1<<52/log(2)) = 6497320848556798, 1023<<52 = 4607182418800017408
 *
 * The fudge factors such that exp y <= exp_fast y:
 * >>> ceiling (2^23 * (1 - (log (log 2) + 1)/log 2))
 * 722019
 * >>> ceiling (2^20 * (1 - (log (log 2) + 1)/log 2))
 * 90253
 * >>> ceiling (2^52 * (1 - (log (log 2) + 1)/log 2))
 * 387630818974388
 *
 * The fudge factor such that exp_fast y <= exp y is uniformly -1
 *
 * TODO: perform exponential doubling for pow based on better_exp_fast instead for better accuracy.
 */

/* Schraudolph's published algorithm extended into the least significant bits to avoid the stair step.
 double long long approximation: round 1<<52/log(2) 6497320848556798,
  mask = 0x3ff0000000000000LL = 4607182418800017408LL
 double approximation: round(1<<20/log(2)) = 1512775, 1023<<20 = 1072693248
*/

/* 4607182418800017408 - 387630818974388 = 4606794787981043020

Exponent mask adapted to full 64 bit precision:
>>> 1023 * 2^52
4607182418800017408

The fudge factor for conservative lower bound adapted to full 64 bit precision:
>>> round (2^52 * (1 - (log (log 2) + 1)/log 2))
387630818974388

As a lower bound this is suitable for use when generating Mass and Precision estimates.
*/
double exp_fast_lb(double a) {
  union { double d; long long x; } u;
  u.x = (long long)(6497320848556798LL * a + 4606794787981043020);
  return u.d;
}

/* 4607182418800017408 + 1 */
double exp_fast_ub(double a) {
  union { double d; long long x; } u;
  u.x = (long long)(6497320848556798LL * a + 4607182418800017409);
  return u.d;
}

double exp_fast(double a) {
  union { double d; long long x; } u;
  u.x = (long long)(6497320848556798LL * a + 0x3fef127e83d16f12LL);
  return u.d;
}

double better_exp_fast(double a) {
  union { double d; long long x; } u, v;
  u.x = (long long)(3248660424278399LL * a + 0x3fdf127e83d16f12LL);
  v.x = (long long)(0x3fdf127e83d16f12LL - 3248660424278399LL * a);
  return u.d / v.d;
}

/* Schraudolph's published algorithm */
double exp_fast_schraudolph(double a) {
  union { double d; int x[2]; } u;
  u.x[1] = (int) (1512775 * a + 1072632447);
  u.x[0] = 0;
  return u.d;
}

/* 1065353216 + 1 */
float expf_fast_ub(float a) {
  union { float f; int x; } u;
  u.x = (int) (12102203 * a + 1065353217);
  return u.f;
}

/* Schraudolph's published algorithm with John's constants */
/* 1065353216 - 486411 = 1064866805 */
float expf_fast(float a) {
  union { float f; int x; } u;
  u.x = (int) (12102203 * a + 1064866805);
  return u.f;
}

//  1056478197 
double better_expf_fast(float a) {
  union { float f; int x; } u, v;
  u.x = (long long)(6051102 * a + 1056478197);
  v.x = (long long)(1056478197 - 6051102 * a);
  return u.f / v.f;
}

/* 1065353216 - 722019 */
float expf_fast_lb(float a) {
  union { float f; int x; } u;
  u.x = (int) (12102203 * a + 1064631197);
  return u.f;
}

/* Ankerl's inversion of Schraudolph's published algorithm, converted to explicit multiplication */
double log_fast_ankerl(double a) {
  union { double d; int x[2]; } u = { a };
  return (u.x[1] - 1072632447) * 6.610368362777016e-7; /* 1 / 1512775.0; */
}

double log_fast_ub(double a) {
  union { double d; long long x; } u = { a };
  return (u.x - 4606794787981043020) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */
}

/* Ankerl's inversion of Schraudolph's published algorithm with my constants */
double log_fast(double a) {
  union { double d; long long x; } u = { a };
  return (u.x - 4606921278410026770) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */
}

double log_fast_lb(double a) {
  union { double d; long long x; } u = { a };
  return (u.x - 4607182418800017409) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */
}


/* 1065353216 - 722019 */
float logf_fast_ub(float a) {
  union { float f; int x; } u = { a };
  return (u.x - 1064631197) * 8.262958405176314e-8f; /* 1 / 12102203.0; */
}

/* Ankerl's adaptation of Schraudolph's published algorithm with John's constants */
/* 1065353216 - 486411 = 1064866805 */
float logf_fast(float a) {
  union { float f; int x; } u = { a };
  return (u.x - 1064866805) * 8.262958405176314e-8f; /* 1 / 12102203.0; */
}

/* 1065353216 + 1 */
float logf_fast_lb(float a) {
  union { float f; int x; } u = { a };
  return (u.x - 1065353217) * 8.262958405176314e-8f; /* 1 / 12102203.0 */
}

/* Ankerl's version of Schraudolph's approximation. */
double pow_fast_ankerl(double a, double b) {
  union { double d; int x[2]; } u = { a };
  u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447);
  u.x[0] = 0;
  return u.d;
}

/*
 These constants are based loosely on the following comment off of Ankerl's blog:

 "I have used the same trick for float, not double, with some slight modification to the constants to suite IEEE754 float format. The first constant for float is 1<<23/log(2) and the second is 127<<23 (for double they are 1<<20/log(2) and 1023<<20)." -- John
*/

/* 1065353216 + 1      = 1065353217 ub */
/* 1065353216 - 486411 = 1064866805 min RMSE */
/* 1065353216 - 722019 = 1064631197 lb */
float powf_fast(float a, float b) {
  union { float d; int x; } u = { a };
  u.x = (int)(b * (u.x - 1064866805) + 1064866805);
  return u.d;
}

float powf_fast_lb(float a, float b) {
  union { float d; int x; } u = { a };
  u.x = (int)(b * (u.x - 1065353217) + 1064631197);
  return u.d;
}

float powf_fast_ub(float a, float b) {
  union { float d; int x; } u = { a };
  u.x = (int)(b * (u.x - 1064631197) + 1065353217);
  return u.d;
}

/*
  Now that 64 bit arithmetic is cheap we can (try to) improve on Ankerl's algorithm.

 double long long approximation: round 1<<52/log(2) 6497320848556798,
  mask = 0x3ff0000000000000LL = 4607182418800017408LL

>>> round (2**52 * log (3 / (8 * log 2) + 1/2) / log 2 - 1/2)
261140389990638
>>> 0x3ff0000000000000 - round (2**52 * log (3 / (8 * log 2) + 1/2) / log 2 - 1/2)
4606921278410026770

*/

double pow_fast_ub(double a, double b) {
  union { double d; long long x; } u = { a };
  u.x = (long long)(b * (u.x - 4606794787981043020LL) + 4607182418800017409LL);
  return u.d;
}

double pow_fast(double a, double b) {
  union { double d; long long x; } u = { a };
  u.x = (long long)(b * (u.x - 4606921278410026770LL) + 4606921278410026770LL);
  return u.d;
}

double pow_fast_lb(double a, double b) {
  union { double d; long long x; } u = { a };
  u.x = (long long)(b * (u.x - 4607182418800017409LL) + 4606794787981043020LL);
  return u.d;
}

/* should be much more precise with large b, still ~3.3x faster. */
double pow_fast_precise_ankerl(double a, double b) {
  int flipped = 0;
  if (b < 0) {
    flipped = 1;
    b = -b;
  }

  /* calculate approximation with fraction of the exponent */
  int e = (int) b;
  union { double d; int x[2]; } u = { a };
  u.x[1] = (int)((b - e) * (u.x[1] - 1072632447) + 1072632447);
  u.x[0] = 0;

  double r = 1.0;
  while (e) {
    if (e & 1) {
      r *= a;
    }
    a *= a;
    e >>= 1;
  }

  r *= u.d;
  return flipped ? 1.0/r : r;
}

/* should be much more precise with large b, still ~3.3x faster. */
double pow_fast_precise(double a, double b) {
  int flipped = 0;
  if (b < 0) {
    flipped = 1;
    b = -b;
  }

  /* calculate approximation with fraction of the exponent */
  int e = (int) b;
  double d = exp_fast(b - e);

  double r = 1.0;
  while (e) {
    if (e & 1) r *= a;
    a *= a;
    e >>= 1;
  }

  r *= d;
  return flipped ? 1.0/r : r;
}

double better_pow_fast_precise(double a, double b) {
  int flipped = 0;
  if (b < 0) {
    flipped = 1;
    b = -b;
  }

  /* calculate approximation with fraction of the exponent */
  int e = (int) b;
  double d = better_exp_fast(b - e);

  double r = 1.0;
  while (e) {
    if (e & 1) r *= a;
    a *= a;
    e >>= 1;
  }

  r *= d;
  return flipped ? 1.0/r : r;
}


/* should be much more precise with large b */
float powf_fast_precise(float a, float b) {
  int flipped = 0;
  if (b < 0) {
    flipped = 1;
    b = -b;
  }

  /* calculate approximation with fraction of the exponent */
  int e = (int) b;
  union { float f; int x; } u = { a };
  u.x = (int)((b - e) * (u.x - 1065353216) + 1065353216);

  float r = 1.0f;
  while (e) {
    if (e & 1) {
      r *= a;
    }
    a *= a;
    e >>= 1;
  }

  r *= u.f;
  return flipped ? 1.0f/r : r;
}

/* should be much more precise with large b */
float better_powf_fast_precise(float a, float b) {
  int flipped = 0;
  if (b < 0) {
    flipped = 1;
    b = -b;
  }

  /* calculate approximation with fraction of the exponent */
  int e = (int) b;
  float f = better_expf_fast(b - e);

  float r = 1.0f;
  while (e) {
    if (e & 1) {
      r *= a;
    }
    a *= a;
    e >>= 1;
  }

  r *= f;
  return flipped ? 1.0f/r : r;
}