expm1 function

double expm1(double x, { List<double>? hiPrecOut, })

Compute exp(x) - 1.

Implementation

double expm1(double x, {List<double>? hiPrecOut}) {
  if (x != x || x == 0.0) {
    // NaN or zero
    return x;
  }

  if (x <= -1.0 || x >= 1.0) {
    // If not between +/- 1.0
    //return exp(x) - 1.0;
    final hiPrec = List<double>.filled(2, 0.0);
    expFast(x, 0.0, hiPrec);
    if (x > 0.0) {
      return -1.0 + hiPrec[0] + hiPrec[1];
    } else {
      final ra = -1.0 + hiPrec[0];
      var rb = -(ra + 1.0 - hiPrec[0]);
      rb += hiPrec[1];
      return ra + rb;
    }
  }

  double baseA;
  double baseB;
  double epsilon;
  var negative = false;

  if (x < 0.0) {
    x = -x;
    negative = true;
  }

  {
    var intFrac = (x * 1024.0).toInt();
    var tempA = EXP_FRAC_TABLE_A[intFrac] - 1.0;
    var tempB = EXP_FRAC_TABLE_B[intFrac];

    var temp = tempA + tempB;
    tempB = -(temp - tempA - tempB);
    tempA = temp;

    temp = tempA * HEX_40000000;
    baseA = tempA + temp - temp;
    baseB = tempB + (tempA - baseA);

    epsilon = x - intFrac / 1024.0;
  }

  /// Compute expm1(epsilon)
  var zb = 0.008336750013465571;
  zb = zb * epsilon + 0.041666663879186654;
  zb = zb * epsilon + 0.16666666666745392;
  zb = zb * epsilon + 0.49999999999999994;
  zb *= epsilon;
  zb *= epsilon;

  var za = epsilon;
  var temp = za + zb;
  zb = -(temp - za - zb);
  za = temp;

  temp = za * HEX_40000000;
  temp = za + temp - temp;
  zb += za - temp;
  za = temp;

  /// Combine the parts.   `expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b)`
  var ya = za * baseA;
  //double yb = za*baseB + zb*baseA + zb*baseB;
  temp = ya + za * baseB;
  var yb = -(temp - ya - za * baseB);
  ya = temp;

  temp = ya + zb * baseA;
  yb += -(temp - ya - zb * baseA);
  ya = temp;

  temp = ya + zb * baseB;
  yb += -(temp - ya - zb * baseB);
  ya = temp;

  //ya = ya + za + baseA;
  //yb = yb + zb + baseB;
  temp = ya + baseA;
  yb += -(temp - baseA - ya);
  ya = temp;

  temp = ya + za;
  //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
  yb += -(temp - ya - za);
  ya = temp;

  temp = ya + baseB;
  //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
  yb += -(temp - ya - baseB);
  ya = temp;

  temp = ya + zb;
  //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
  yb += -(temp - ya - zb);
  ya = temp;

  if (negative) {
    /// Compute `expm1(-x) = -expm1(x) / (expm1(x) + 1)`
    var denom = 1.0 + ya;
    var denomr = 1.0 / denom;
    var denomb = -(denom - 1.0 - ya) + yb;
    var ratio = ya * denomr;
    temp = ratio * HEX_40000000;
    final ra = ratio + temp - temp;
    var rb = ratio - ra;

    temp = denom * HEX_40000000;
    za = denom + temp - temp;
    zb = denom - za;

    rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;

    // f(x) = x/1+x
    // Compute f'(x)
    // Product rule:  d(uv) = du*v + u*dv
    // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
    // d(1/x) = -1/(x*x)
    // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
    // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))

    // Adjust for yb
    rb += yb * denomr; // numerator
    rb += -ya * denomb * denomr * denomr; // denominator

    // negate
    ya = -ra;
    yb = -rb;
  }

  if (hiPrecOut != null) {
    hiPrecOut[0] = ya;
    hiPrecOut[1] = yb;
  }

  return ya + yb;
}