gammaLowerRegularized function

double gammaLowerRegularized (double a, double x)

Returns the lower incomplete regularized gamma function P(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0.

Implementation

double gammaLowerRegularized(double a, double x) {
  const double epsilon = 0.000000000000001;
  const double big = 4503599627370496.0;
  const double bigInv = 2.22044604925031308085e-16;

  if (a < 0.0) {
    throw ArgumentError.value(a, 'a', messages.argumentNotNegative);
  }

  if (x < 0.0) {
    throw ArgumentError.value(x, 'x', messages.argumentNotNegative);
  }

  if (almostEqual(a, 0.0)) {
    if (almostEqual(x, 0.0)) {
      //use right hand limit value because so that regularized upper/lower gamma definition holds.
      return 1.0;
    }

    return 1.0;
  }

  if (almostEqual(x, 0.0)) {
    return 0.0;
  }

  double ax = (a * stdmath.log(x)) - x - gammaLn(a);
  if (ax < -709.78271289338399) {
    return a < x ? 1.0 : 0.0;
  }

  if (x <= 1 || x <= a) {
    double r2 = a;
    double c2 = 1.0;
    double ans2 = 1.0;

    do {
      r2 = r2 + 1;
      c2 = c2 * x / r2;
      ans2 += c2;
    } while ((c2 / ans2) > epsilon);

    return stdmath.exp(ax) * ans2 / a;
  }

  int c = 0;
  double y = 1 - a;
  double z = x + y + 1;

  double p3 = 1.0;
  double q3 = x;
  double p2 = x + 1;
  double q2 = z * x;
  double ans = p2 / q2;

  double error;

  do {
    c++;
    y += 1;
    z += 2;
    double yc = y * c;

    double p = (p2 * z) - (p3 * yc);
    double q = (q2 * z) - (q3 * yc);

    if (q != 0) {
      double nextans = p / q;
      error = ((ans - nextans) / nextans).abs();
      ans = nextans;
    } else {
      // zero div, skip
      error = 1.0;
    }

    // shift
    p3 = p2;
    p2 = p;
    q3 = q2;
    q2 = q;

    // normalize fraction when the numerator becomes large
    if (p.abs() > big) {
      p3 *= bigInv;
      p2 *= bigInv;
      q3 *= bigInv;
      q2 *= bigInv;
    }
  } while (error > epsilon);

  return 1.0 - (stdmath.exp(ax) * ans);
}