exp function

double exp(double x)

Optimized Exponential function.

• x is bounded to -87 to 87, since bigger values won't have enough precision when stored at Float32x4.

NOTE: For higher precision use expHighPrecision.

Implementation

double exp(double x) {
  x = x.clamp(-87, 87);

  var xInt = _expXToInt(x);
  final xIdx = expIntTableMaxIndex + xInt;

  // Lookup exp(floor(x)).
  // intPartA will have the upper 22 bits, intPartB will have the lower
  // 52 bits.
  final intPartA = expIntTableA[xIdx];
  final intPartB = expIntTableB[xIdx];

  // Get the fractional part of x, find the greatest multiple of 2^-10 less than
  // x and look up the exp function of it.
  // fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
  final fracInt = ((x - xInt) * 1024.0).toInt();

  final fracPartA = expFracTableA[fracInt];
  final fracPartB = expFracTableB[fracInt];

  // epsilon is the difference in x from the nearest multiple of 2^-10.  It
  // has a value in the range 0 <= epsilon < 2^-10.
  // Do the subtraction from x as the last step to avoid possible
  // loss of precision.
  final epsilon = x - (xInt + fracInt / 1024.0);

  // Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
  // full double precision (52 bits).  Since z < 2^-10, we will have
  // 62 bits of precision when combined with the contant 1.  This will be
  // used in the last addition below to get proper rounding.

  // Remez generated polynomial.  Converges on the interval [0, 2^-10], error
  // is less than 0.5 ULP
  var o = 0.04168701738764507;
  o = o * epsilon + 0.1666666505023083;
  o = o * epsilon + 0.5000000000042687;
  o = o * epsilon + 1.0;
  o = o * epsilon + -3.940510424527919E-20;

  // Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
  // expansion.
  // tempA is exact since intPartA and intPartB only have 22 bits each.
  // tempB will have 52 bits of precision.
  final tempA = intPartA * fracPartA;
  final tempB =
      intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;

  // Compute the result.  (1+z)(tempA+tempB).  Order of operations is
  // important.  For accuracy add by increasing size.  tempA is exact and
  // much larger than the others.  If there are extra bits specified from the
  // pow() function, use them.
  final tempC = tempB + tempA;
  var result = tempC * o + tempB + tempA;

  return result;
}