gammaUpperRegularized function

double gammaUpperRegularized (double a, double x)

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

Implementation

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

  if (x < 1.0 || x <= a) {
    return 1.0 - gammaLowerRegularized(a, x);
  }

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

  ax = stdmath.exp(ax);
  double t;
  double y = 1 - a;
  double z = x + y + 1;
  double c = 0.0;
  double pkm2 = 1.0;
  double qkm2 = x;
  double pkm1 = x + 1;
  double qkm1 = z * x;
  double ans = pkm1 / qkm1;
  do {
    c = c + 1;
    y = y + 1;
    z = z + 2;
    double yc = y * c;
    double pk = pkm1 * z - pkm2 * yc;
    double qk = qkm1 * z - qkm2 * yc;
    if (qk != 0) {
      double r = pk / qk;
      t = ((ans - r) / r).abs();
      ans = r;
    } else {
      t = 1.0;
    }

    pkm2 = pkm1;
    pkm1 = pk;
    qkm2 = qkm1;
    qkm1 = qk;

    if (pk.abs() > big) {
      pkm2 = pkm2 * bigInv;
      pkm1 = pkm1 * bigInv;
      qkm2 = qkm2 * bigInv;
      qkm1 = qkm1 * bigInv;
    }
  } while (t > epsilon);

  return ans * ax;
}