sinh function

double sinh(double x)

Compute the hyperbolic sine of a number.

Implementation

double sinh(double x) {
  if (x != x) {
    return x;
  }

  var negate = false;

  // sinh[z] = (exp(z) - exp(-z) / 2

  // for values of z larger than about 20,
  // exp(-z) can be ignored in comparison with exp(z)

  if (x > 20) {
    if (x >= logMaxValue) {
      // Avoid overflow (MATH-905).
      final t = exp(0.5 * x);
      return (0.5 * t) * t;
    } else {
      return 0.5 * exp(x);
    }
  } else if (x < -20) {
    if (x <= -logMaxValue) {
      // Avoid overflow (MATH-905).
      final t = exp(-0.5 * x);
      return (-0.5 * t) * t;
    } else {
      return -0.5 * exp(-x);
    }
  }

  if (x == 0) {
    return x;
  }

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

  double result;

  if (x > 0.25) {
    final hiPrec = List.filled(2, 0.0);
    exp(x, 0.0, hiPrec);

    var ya = hiPrec[0] + hiPrec[1];
    var yb = -(ya - hiPrec[0] - hiPrec[1]);

    var temp = ya * hex40000000;
    final yaa = ya + temp - temp;
    final yab = ya - yaa;

    // recip = 1/y
    final recip = 1.0 / ya;
    temp = recip * hex40000000;
    var recipa = recip + temp - temp;
    var recipb = recip - recipa;

    // Correct for rounding in division
    recipb +=
        (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) *
            recip;
    // Account for yb
    recipb += -yb * recip * recip;

    recipa = -recipa;
    recipb = -recipb;

    // y = y + 1/y
    temp = ya + recipa;
    yb += -(temp - ya - recipa);
    ya = temp;
    temp = ya + recipb;
    yb += -(temp - ya - recipb);
    ya = temp;

    result = ya + yb;
    result *= 0.5;
  } else {
    final hiPrec = List.filled(2, 0.0);
    expm1(x, hiPrec);

    var ya = hiPrec[0] + hiPrec[1];
    var yb = -(ya - hiPrec[0] - hiPrec[1]);

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

    temp = denom * hex40000000;
    final za = denom + temp - temp;
    final zb = denom - za;

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

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

    // y = y - 1/y
    temp = ya + ra;
    yb += -(temp - ya - ra);
    ya = temp;
    temp = ya + rb;
    yb += -(temp - ya - rb);
    ya = temp;

    result = ya + yb;
    result *= 0.5;
  }

  if (negate) {
    result = -result;
  }

  return result;
}