gamma function

double gamma(double x)

Implementation

double gamma(double x) {
  const double _euler =
      0.57721566490153286060651209008240243104215933593992; // A001620

  // Special cases
  if (x.isNegative && Double.fractionPart(x) == 0) {
    return double.nan;
  }

  if ((x.isInfinite && x.isNegative) || x.isNaN) {
    return double.nan;
  }

  if (x.isInfinite) return double.infinity;

  if (x == 0) {
    if (x.isNegative) return double.negativeInfinity;
    return double.infinity;
  }

  double q = x.abs();
  double p = q.floorToDouble();

  if (q > 33) {
    if (x >= 0) {
      List<double> y = _stirling(x);
      return y.first * y.last;
    }
    // Note: x is negative but (checked above) not a negative integer,
    // so x must be small enough to be in range for conversion to int64.
    // If |x| were >= 2⁶³ it would have to be an integer.
    int signgam = 1;
    {
      final int ip = p.toInt();
      if (ip & 1 == 0) {
        signgam = -1;
      }
    }
    double z = q - p;
    if (z > 0.5) {
      p = p + 1;
      z = q - p;
    }
    z = q * math.sin(math.pi * z);
    if (z == 0) {
      return signgam > 0 ? double.infinity : double.negativeInfinity;
    }
    final List<double> sq = _stirling(q);
    final double absz = z.abs();
    final double d = absz * sq.first * sq.last;
    if (d.isInfinite) {
      z = math.pi / absz / sq.first / sq.last;
    } else {
      z = math.pi / d;
    }
    return signgam * z;
  }

  // Reduce argument
  double z = 1.0;
  while (x >= 3) {
    x = x - 1;
    z = z * x;
  }
  while (x < 0) {
    if (x > -1e-09) {
      if (x == 0) {
        return double.infinity;
      }
      return z / ((1 + _euler * x) * x);
    }
    z = z / x;
    x = x + 1;
  }
  while (x < 2) {
    if (x < 1e-09) {
      if (x == 0) {
        return double.infinity;
      }
      return z / ((1 + _euler * x) * x);
    }
    z = z / x;
    x = x + 1;
  }

  if (x == 2) {
    return z;
  }

  x = x - 2;
  p = (((((x * _gamP[0] + _gamP[1]) * x + _gamP[2]) * x + _gamP[3]) * x +
                      _gamP[4]) *
                  x +
              _gamP[5]) *
          x +
      _gamP[6];
  q = ((((((x * _gamQ[0] + _gamQ[1]) * x + _gamQ[2]) * x + _gamQ[3]) * x +
                              _gamQ[4]) *
                          x +
                      _gamQ[5]) *
                  x +
              _gamQ[6]) *
          x +
      _gamQ[7];
  return z * p / q;
}