Regularized Incomplete Gamma Function P(a,x) = <i><big>&#8747;</big><sub><small>0</small></sub><sup><small>x</small></sup> e<sup>-t</sup> t<sup>(a-1)</sup> dt</i>.
Series representation of the function - valid for x < a + 1
322023
private static double regularizedIncompleteGammaSeries(double a, double x) {
        if (a < 0.0 || x < 0.0 || x >= a + 1) {
            throw new IllegalArgumentException(String.format("Invalid a = %f, x = %f", a, x));
        }

        int i = 0;
        double igf = 0.0;
        boolean check = true;

        double acopy = a;
        double sum = 1.0 / a;
        double incr = sum;
        double loggamma = lgamma(a);

        while (check) {
            ++i;
            ++a;
            incr *= x / a;
            sum += incr;
            if (Math.abs(incr) < Math.abs(sum) * INCOMPLETE_GAMMA_EPSILON) {
                igf = sum * Math.exp(-x + acopy * Math.log(x) - loggamma);
                check = false;
            }
            if (i >= INCOMPLETE_GAMMA_MAX_ITERATIONS) {
                check = false;
                igf = sum * Math.exp(-x + acopy * Math.log(x) - loggamma);
                logger.error("Gamma.regularizedIncompleteGammaSeries: Maximum number of iterations wes exceeded");
            }
        }
        return igf;
    }
765130
private function gamma($x)
    {
        $gamma = 0.577215664901532860606512090; # Euler's gamma constant
        if ($x < 0.001) {
            return 1.0/($x*(1.0 + $gamma*$x));
        }
        if ($x < 12.0) {
            # The algorithm directly approximates gamma over (1,2) and uses
            # reduction identities to reduce other arguments to this interval.
            $y = $x;
            $n = 0;
            $arg_was_less_than_one = ($y < 1.0);
            # Add or subtract integers as necessary to bring y into (1,2)
            # Will correct for this below
            if ($arg_was_less_than_one) {
                $y += 1.0;
            } else {
                $n = floor($y) - 1;  # will use n later
                $y -= $n;
            }
            # numerator coefficients for approximation over the interval (1,2)
            $p =
            array(
                -1.71618513886549492533811E+0,
                 2.47656508055759199108314E+1,
                -3.79804256470945635097577E+2,
                 6.29331155312818442661052E+2,
                 8.66966202790413211295064E+2,
                -3.14512729688483675254357E+4,
                -3.61444134186911729807069E+4,
                 6.64561438202405440627855E+4
            );

            # denominator coefficients for approximation over the interval (1,2)
            $q =
            array(
                -3.08402300119738975254353E+1,
                 3.15350626979604161529144E+2,
                -1.01515636749021914166146E+3,
                -3.10777167157231109440444E+3,
                 2.25381184209801510330112E+4,
                 4.75584627752788110767815E+3,
                -1.34659959864969306392456E+5,
                -1.15132259675553483497211E+5
            );

            $num = 0.0;
            $den = 1.0;

            $z = $y - 1;
            for ($i = 0; $i < 8; $i++) {
                $num = ($num + $p[$i])*$z;
                $den = $den*$z + $q[$i];
            }
            $result = $num/$den + 1.0;

            # Apply correction if argument was not initially in (1,2)
            if ($arg_was_less_than_one) {
                # Use identity gamma(z) = gamma(z+1)/z
                # The variable "result" now holds gamma of the original y + 1
                # Thus we use y-1 to get back the orginal y.
                $result /= ($y-1.0);
            } else {
                # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
                for ($i = 0; $i < $n; $i++) {
                    $result *= $y++;
                }
            }

            return $result;
        }

        ###########################################################################
        # Third interval: [12, infinity)

        if ($x > 171.624) {
            # Correct answer too large to display.

            return Double.POSITIVE_INFINITY;
        }

        return exp($this->log_gamma($x));
    }
1082636
def get_variable_errors(X, L, LX, Z, U, step_g, e_rel, e_abs=0):
    
    n = X.size
    p = Z.size
    e_pri2 = np.sqrt(p)*e_abs/L.spectral_norm + e_rel*np.max([l2(LX), l2(Z)])
    if step_g is not None:
        e_dual2 = np.sqrt(n)*e_abs/L.spectral_norm + e_rel*l2(L.T.dot(U)/step_g)
    else:
        e_dual2 = np.sqrt(n)*e_abs/L.spectral_norm + e_rel*l2(L.T.dot(U))
    return e_pri2, e_dual2
150076
function findx0SmallA(a, p, q) {
      var B, u, temp;

      B = q * gamma(a);
      if (B > 0.6 || B >= 0.45 && a >= 0.3) {
         // use 21
         u = B * q > 1e-8 ? Math.pow(p * gamma(a + 1), 1 / a)
                          : Math.exp(-q / a - eulerGamma);
         return u / (1 - u / (a + 1));
      }
      if (a < 0.3 && B >= 0.35 && B <= 0.6) {
         // use 22
         temp = Math.exp(-eulerGamma - B);
         u = temp * Math.exp(temp);
         return temp * Math.exp(u);
      }
      temp = -Math.log(B);
      u = temp - (1 - a) * Math.log(temp);
      if (B >= 0.15) {
         // use 23, with temp for y and u for v
         return temp - (1 - a) * Math.log(u) -
                  Math.log(1 + (1 - a) / (1 + u));
      }
      if (B > 0.01) {
         // use 24, with temp for y and u for v
         return temp - (1 - a) * Math.log(u) - Math.log(
            (u * u + 2 * (3 - a) * u + (2 - a) * (3 - a)) /
            (u * u + (5 - a) * u + 2));
      }
      // use 25, where c1 + ... + c5/y^4 is treated as a polynomial in y^-1
      // using u for c1, temp for y
      return findx0TinyB(a, temp);
   }
848550
def GaP_Eg_Gamma(self, **kwargs):
    '''
    
    '''
    T = kwargs.get('T', 300.)
    if T < 1e-4:
        return self.Eg_Gamma_0()
    return self.Eg_Gamma_0() + 0.1081 * (1 - 1. / tanh(164. / T))  # eV
600915
function brcomp(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      if (Math.min(a, b) < 8) {
         return a * Math.log(x) + b * Math.log(y) - lbeta(a, b);
      }
      return 0.5 * (Math.log(a) + Math.log(b) - C.log2pi - Math.log(a + b)) -
             a * phi(1 - lambda / a) - b * phi(1 + lambda / b) -
             bcorr(a, b);
   }
790101
function bpser(a, b, x) {
      var c, ser;

      c = 1;
      ser = a * series(function(i) {
         if (i === 0) { return 0; }
         c = c * (1 - b / i) * x;
         return c / (a + i);
      });
      ser = ser > -1 ? log1p(ser) : -Infinity;
      return -lbeta(a, b) + a * Math.log(x) - Math.log(a) + ser;
   }
1052097
function p_t(df, t) {
        var c2, s, p, i;

        c2 = df / (df + t * t);
        s = Math.sqrt(1.0 - c2);

        if (t < 0.0) { s = -s; }

        p = 0.0;
        i = df % 2 + 2;
        while (i <= df) {
            p = p + s;
            s = s * ((i - 1) * c2 / i);
            i = i + 2;
        }

        if ((df & 1) !== 0) {
            return 0.5 + (p * Math.sqrt(c2) + Math.atan(t / Math.sqrt(df))) / Math.PI;
        }

        return (1.0 + p) / 2.0;
    }
671369
function gam1(x) {
      if (x <= -1) { return NaN; }
      if (x < -0.5 || x > 1.5) { return 1 / gamma(x + 1) - 1; }
      if (x <= 0) { return x * (1 + w.evalAt(x)); }
      if (x <= 0.5) { return x * w1.evalAt(x); }
      if (x <= 1) { return (x - 1) / x * w.evalAt(x - 1); }
      return (x - 1) / x * (w1.evalAt(x - 1) - 1);
   }
352889
function bup(a, b, x, y, n) {
      var logSeries;

      logSeries = Math.log(series(function(i, v) {
            if (i === 0) { return 1; }
            return v * x * (1 + (b - 1) / (a + i));
         }, n));
      return a * Math.log(x) + b * Math.log(y) - lbeta(a, b) - Math.log(a) +
         logSeries;
   }
292189
function gamma1_64(x)
	{
		var rotr19 = rotr_64(x, 19), rotr61 = rotr_64(x, 61),
			shr6 = shr_64(x, 6);

		return new Int_64(
				rotr19.highOrder ^ rotr61.highOrder ^ shr6.highOrder,
				rotr19.lowOrder ^ rotr61.lowOrder ^ shr6.lowOrder
			);
	}
382572
function ecp_proj(R, N, p, l) {
	var x = N[0], y = N[1], z = N[2];
	var xr = R[0], yr = R[1], zr = R[2];
	var t0 = _ecp_t0, t1 = _ecp_t1, t2 = _ecp_t2;
	fp_mul(t2, t0, fp_sqr(t1, fp_inv(t0, z, p, l), p, l), p, l);
	fp_mul(xr, x, t1, p, l);
	fp_mul(yr, y, t2, p, l);
	mpn_zero(zr, l), zr[0] = 1;
	return R;
}
1141866
def gammalnStirling(z):
    
    return (0.5 * (np.log(2. * np.pi) - np.log(z))) \
           + (z * (np.log(z + (1. / ((12. * z) - (1. / (10. * z))))) - 1.))
574169
private double[][] estimateGamma(double[][][] xi) {
        double[][] gamma = new double[xi.length + 1][numStates];

        for (int t = 0; t < xi.length; t++) {
            for (int i = 0; i < numStates; i++) {
                for (int j = 0; j < numStates; j++) {
                    gamma[t][i] += xi[t][i][j];
                }
            }
        }

        int n = xi.length - 1;
        for (int j = 0; j < numStates; j++) {
            for (int i = 0; i < numStates; i++) {
                gamma[xi.length][j] += xi[n][i][j];
            }
        }

        return gamma;
    }
88943
function partial( lambda ) {

	/**
	* FUNCTION: cdf( x )
	*	Evaluates the cumulative distribution function (CDF) for a Poisson distribution.
	*
	* @private
	* @param {Number} x - input value
	* @returns {Number} evaluated CDF
	*/
	return function cdf( x ) {
		if ( x < 0 ) {
			return 0;
		}
		return gammainc( lambda, floor( x ) + 1, {
			'tail': 'upper',
			'regularized': true
		});
	};
}
867563
function loggamma(x) {
        var v = 1.0,
            w = 1.0 / (x * x),
            ret;
        while (x < N) {
            v = v * x;
            x = x + 1;
        }
        ret = B16 / (16 * 15);
        ret = ret * w + B14 / (14 * 13);
        ret = ret * w + B12 / (12 * 11);
        ret = ret * w + B10 / (10 *  9);
        ret = ret * w + B8  / (8 *  7);
        ret = ret * w + B6  / (6 *  5);
        ret = ret * w + B4  / (4 *  3);
        ret = ret * w + B2  / (2 *  1);
        ret = ret / x + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);
        return ret;
    }
1370914
def lpgamma(p, a):
    
    res = p * (p - 1.0) / 4.0 * log(pi)
    for ii in range(1, p+1):
        res += gammaln(a + (1.0 - ii) / 2.0)

    return res