expHighPrecision function

double expHighPrecision(double x, [ double extra = 0.0, List<double>? highPrecision ])

Internal helper method for exponential function.

x is the original argument of the exponential function. extra bits of precision on input (To Be Confirmed). highPrecision extra bits of precision on output (To Be Confirmed)

Implementation

double expHighPrecision(
  double x, [
  double extra = 0.0,
  List<double>? highPrecision,
]) {
  double intPartA;
  double intPartB;
  int intVal;

  if (x != x || x.isInfinite) {
    // `exp(NaN) = NaN`, `exp(+inf) = +inf` and `exp(-inf) = 0`:
    var result = x != x ? double.nan : (x > 0 ? double.infinity : 0.0);

    if (highPrecision != null) {
      highPrecision[0] = result;
      highPrecision[1] = 0.0;
    }

    return result;
  }

  // Lookup exp(floor(x)).
  // intPartA will have the upper 22 bits, intPartB will have the lower
  // 52 bits.
  if (x < 0.0) {
    intVal = -x.toInt();

    if (intVal > 746) {
      if (highPrecision != null) {
        highPrecision[0] = 0.0;
        highPrecision[1] = 0.0;
      }
      return 0.0;
    }

    if (intVal > 709) {
      // This will produce a subnormal output
      final result =
          expHighPrecision(x + 40.19140625, extra, highPrecision) /
          285040095144011776.0;
      if (highPrecision != null) {
        highPrecision[0] /= 285040095144011776.0;
        highPrecision[1] /= 285040095144011776.0;
      }
      return result;
    }

    if (intVal == 709) {
      // exp(1.494140625) is nearly a machine number...
      final result =
          expHighPrecision(x + 1.494140625, extra, highPrecision) /
          4.455505956692756620;
      if (highPrecision != null) {
        highPrecision[0] /= 4.455505956692756620;
        highPrecision[1] /= 4.455505956692756620;
      }
      return result;
    }

    intVal++;

    intPartA = expIntTableA[expIntTableMaxIndex - intVal];
    intPartB = expIntTableB[expIntTableMaxIndex - intVal];

    intVal = -intVal;
  } else {
    // NOTE: `x` is finite here, the infinities are handled at the top.
    intVal = x.toInt();

    if (intVal > 709) {
      if (highPrecision != null) {
        highPrecision[0] = double.infinity;
        highPrecision[1] = 0.0;
      }
      return double.infinity;
    }

    intPartA = expIntTableA[expIntTableMaxIndex + intVal];
    intPartB = expIntTableB[expIntTableMaxIndex + intVal];
  }

  // 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 intFrac = ((x - intVal) * 1024.0).toInt();
  final fracPartA = expFracTableA[intFrac];
  final fracPartB = expFracTableB[intFrac];

  // 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 percison.
  final epsilon = x - (intVal + intFrac / 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 z = 0.04168701738764507;
  z = z * epsilon + 0.1666666505023083;
  z = z * epsilon + 0.5000000000042687;
  z = z * epsilon + 1.0;
  z = z * 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;
  double result;
  if (extra != 0.0) {
    result = tempC * extra * z + tempC * extra + tempC * z + tempB + tempA;
  } else {
    result = tempC * z + tempB + tempA;
  }

  if (highPrecision != null) {
    // If requesting high precision
    highPrecision[0] = tempA;
    highPrecision[1] = tempC * extra * z + tempC * extra + tempC * z + tempB;
  }

  return result;
}