# tablepress/3.4/libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php

TablePress – Tables in WordPress made easy, version 3.4. 491 lines.

- Page: https://pluginprobe.com/plugins/tablepress/3.4/code/libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php
- Raw: https://pluginprobe.com/plugins/tablepress/3.4/raw/libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php
- Modified: 2026-09-30T03:50:46+00:00

Line numbers below start at 1. Link to a line or a range by appending a fragment to the
page URL, for example `https://pluginprobe.com/plugins/tablepress/3.4/code/libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php#L10-L20`.

```php
<?php

namespace TablePress\PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;

use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Functions;
use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;

abstract class GammaBase
{
	private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;

	private const EPS = 2.22e-16;

	protected const MAX_VALUE = 1.2e308;

	private const SQRT2PI = 2.5066282746310005024157652848110452530069867406099;

	private const MAX_ITERATIONS = 256;

	protected static function calculateDistribution(float $value, float $a, float $b, bool $cumulative): float
	{
		if ($cumulative) {
			return self::regularizedGammaP($a, $value / $b);
		}

		if ($value == 0.0) {
			if ($a == 1.0) {
				return 1.0 / $b;
			}

			return ($a < 1.0) ? INF : 0.0;
		}

		// Log domain, so large shape parameters cannot overflow Gamma(a).
		return exp(($a - 1.0) * log($value) - $value / $b - $a * log($b) - self::logGamma($a));
	}

	/** @return float|string */
	protected static function calculateInverse(float $probability, float $alpha, float $beta)
	{
		$xLo = 0;
		$xHi = $alpha * $beta * 5;

		// Extend the upper bound while it does not yet bracket the root, so a
		// tail quantile beyond alpha*beta*5 is no longer clamped to it. Stop if
		// the CDF stops increasing (series approximation past its usable range)
		// and keep the original bound rather than expanding into that region.
		$xHiBase = $xHi;
		$cdfHi = self::calculateDistribution($xHi, $alpha, $beta, true);
		while ($cdfHi < $probability) {
			$xHiNext = $xHi * 2;
			$cdfNext = self::calculateDistribution($xHiNext, $alpha, $beta, true);
			if ($cdfNext <= $cdfHi) {
				$xHi = $xHiBase;

				break;
			}
			$xHi = $xHiNext;
			$cdfHi = $cdfNext;
		}

		$dx = 1024;
		$x = $xNew = 1;
		$i = 0;

		while ((abs($dx) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
			// Apply Newton-Raphson step
			$result = self::calculateDistribution($x, $alpha, $beta, true);
			$error = $result - $probability;

			if ($error == 0.0) {
				$dx = 0;
			} elseif ($error < 0.0) {
				$xLo = $x;
			} else {
				$xHi = $x;
			}

			$pdf = self::calculateDistribution($x, $alpha, $beta, false);
			// Avoid division by zero
			if ($pdf !== 0.0) {
				$dx = $error / $pdf;
				$xNew = $x - $dx;
			}

			// If the NR fails to converge (which for example may be the
			// case if the initial guess is too rough) we apply a bisection
			// step to determine a more narrow interval around the root.
			if (($xNew < $xLo) || ($xNew > $xHi) || ($pdf == 0.0)) {
				$xNew = ($xLo + $xHi) / 2;
				$dx = $xNew - $x;
			}
			$x = $xNew;
		}

		if ($i === self::MAX_ITERATIONS) {
			return ExcelError::NA();
		}

		return $x;
	}

	//
	//    Implementation of the incomplete Gamma function
	//
	public static function incompleteGamma(float $a, float $x): float
	{
		// Unregularized lower incomplete gamma; kept for backward compatibility.
		return self::regularizedGammaP($a, $x) * self::gammaValue($a);
	}

	/**
	 * Regularized lower incomplete gamma P(a,x) = gamma(a,x) / Gamma(a).
	 * Series for x < a+1, else the complement of the continued fraction.
	 */
	public static function regularizedGammaP(float $a, float $x): float
	{
		if ($x <= 0.0 || $a <= 0.0) {
			return 0.0;
		}
		if ($x < $a + 1.0) {
			return self::gammaSeries($a, $x);
		}

		return 1.0 - self::gammaContinuedFraction($a, $x);
	}

	/**
	 * Regularized upper incomplete gamma Q(a,x) = 1 - P(a,x).
	 * Continued fraction for x >= a+1 keeps the right tail free of cancellation.
	 */
	public static function regularizedGammaQ(float $a, float $x): float
	{
		if ($x <= 0.0 || $a <= 0.0) {
			return 1.0;
		}
		if ($x < $a + 1.0) {
			return 1.0 - self::gammaSeries($a, $x);
		}

		return self::gammaContinuedFraction($a, $x);
	}

	// Near x ~ a both expansions need O(sqrt(a)) terms to reach EPS.
	private static function incompleteGammaIterations(float $a): int
	{
		return max(self::MAX_ITERATIONS, (int) ceil(10.0 * sqrt($a)));
	}

	// P(a,x) by its series representation (Numerical Recipes gser).
	private static function gammaSeries(float $a, float $x): float
	{
		$maxIterations = self::incompleteGammaIterations($a);
		$gln = self::logGamma($a);
		$ap = $a;
		$sum = 1.0 / $a;
		$del = $sum;
		for ($i = 1; $i <= $maxIterations; ++$i) {
			++$ap;
			$del *= $x / $ap;
			$sum += $del;
			if (abs($del) < abs($sum) * self::EPS) {
				break;
			}
		}

		return $sum * exp(-$x + $a * log($x) - $gln);
	}

	// Q(a,x) by its continued fraction representation (Numerical Recipes gcf).
	private static function gammaContinuedFraction(float $a, float $x): float
	{
		$maxIterations = self::incompleteGammaIterations($a);
		$fpMin = 1.0e-300;
		$gln = self::logGamma($a);
		$b = $x + 1.0 - $a;
		$c = 1.0 / $fpMin;
		$d = 1.0 / $b;
		$h = $d;
		for ($i = 1; $i <= $maxIterations; ++$i) {
			$an = -$i * ($i - $a);
			$b += 2.0;
			$d = $an * $d + $b;
			if (abs($d) < $fpMin) {
				$d = $fpMin;
			}
			$c = $b + $an / $c;
			if (abs($c) < $fpMin) {
				$c = $fpMin;
			}
			$d = 1.0 / $d;
			$del = $d * $c;
			$h *= $del;
			if (abs($del - 1.0) < self::EPS) {
				break;
			}
		}

		return $h * exp(-$x + $a * log($x) - $gln);
	}

	private const GAMMA_VALUE_P0 = 1.000000000190015;
	private const GAMMA_VALUE_P = [
		1 => 76.18009172947146,
		2 => -86.50532032941677,
		3 => 24.01409824083091,
		4 => -1.231739572450155,
		5 => 1.208650973866179e-3,
		6 => -5.395239384953e-6,
	];

	//
	//    Implementation of the Gamma function
	//
	public static function gammaValue(float $value): float
	{
		if ($value == 0.0) {
			return 0;
		}

		$y = $x = $value;
		$tmp = $x + 5.5;
		$tmp -= ($x + 0.5) * log($tmp);

		$summer = self::GAMMA_VALUE_P0;
		for ($j = 1; $j <= 6; ++$j) {
			$summer += (self::GAMMA_VALUE_P[$j] / ++$y);
		}

		return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x));
	}

	private const LG_D1 = -0.5772156649015328605195174;

	private const LG_D2 = 0.4227843350984671393993777;

	private const LG_D4 = 1.791759469228055000094023;

	private const LG_P1 = [
		4.945235359296727046734888,
		201.8112620856775083915565,
		2290.838373831346393026739,
		11319.67205903380828685045,
		28557.24635671635335736389,
		38484.96228443793359990269,
		26377.48787624195437963534,
		7225.813979700288197698961,
	];

	private const LG_P2 = [
		4.974607845568932035012064,
		542.4138599891070494101986,
		15506.93864978364947665077,
		184793.2904445632425417223,
		1088204.76946882876749847,
		3338152.967987029735917223,
		5106661.678927352456275255,
		3074109.054850539556250927,
	];

	private const LG_P4 = [
		14745.02166059939948905062,
		2426813.369486704502836312,
		121475557.4045093227939592,
		2663432449.630976949898078,
		29403789566.34553899906876,
		170266573776.5398868392998,
		492612579337.743088758812,
		560625185622.3951465078242,
	];

	private const LG_Q1 = [
		67.48212550303777196073036,
		1113.332393857199323513008,
		7738.757056935398733233834,
		27639.87074403340708898585,
		54993.10206226157329794414,
		61611.22180066002127833352,
		36351.27591501940507276287,
		8785.536302431013170870835,
	];

	private const LG_Q2 = [
		183.0328399370592604055942,
		7765.049321445005871323047,
		133190.3827966074194402448,
		1136705.821321969608938755,
		5267964.117437946917577538,
		13467014.54311101692290052,
		17827365.30353274213975932,
		9533095.591844353613395747,
	];

	private const LG_Q4 = [
		2690.530175870899333379843,
		639388.5654300092398984238,
		41355999.30241388052042842,
		1120872109.61614794137657,
		14886137286.78813811542398,
		101680358627.2438228077304,
		341747634550.7377132798597,
		446315818741.9713286462081,
	];

	private const LG_C = [
		-0.001910444077728,
		8.4171387781295e-4,
		-5.952379913043012e-4,
		7.93650793500350248e-4,
		-0.002777777777777681622553,
		0.08333333333333333331554247,
		0.0057083835261,
	];

	// Rough estimate of the fourth root of logGamma_xBig
	private const LG_FRTBIG = 2.25e76;

	private const PNT68 = 0.6796875;

	// Function cache for logGamma

	private static float $logGammaCacheResult = 0.0;

	private static float $logGammaCacheX = 0.0;

	/**
	 * logGamma function.
	 *
	 * Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
	 *
	 * The natural logarithm of the gamma function. <br />
	 * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz <br />
	 * Applied Mathematics Division <br />
	 * Argonne National Laboratory <br />
	 * Argonne, IL 60439 <br />
	 * <p>
	 * References:
	 * <ol>
	 * <li>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural
	 *     Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.</li>
	 * <li>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.</li>
	 * <li>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.</li>
	 * </ol>
	 * </p>
	 * <p>
	 * From the original documentation:
	 * </p>
	 * <p>
	 * This routine calculates the LOG(GAMMA) function for a positive real argument X.
	 * Computation is based on an algorithm outlined in references 1 and 2.
	 * The program uses rational functions that theoretically approximate LOG(GAMMA)
	 * to at least 18 significant decimal digits. The approximation for X > 12 is from
	 * reference 3, while approximations for X < 12.0 are similar to those in reference
	 * 1, but are unpublished. The accuracy achieved depends on the arithmetic system,
	 * the compiler, the intrinsic functions, and proper selection of the
	 * machine-dependent constants.
	 * </p>
	 * <p>
	 * Error returns: <br />
	 * The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
	 * The computation is believed to be free of underflow and overflow.
	 * </p>
	 *
	 * @version 1.1
	 *
	 * @author Jaco van Kooten
	 *
	 * @return float MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
	 */
	public static function logGamma(float $x): float
	{
		if ($x == self::$logGammaCacheX) {
			return self::$logGammaCacheResult;
		}

		$y = $x;
		if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) {
			if ($y <= self::EPS) {
				$res = -log($y);
			} elseif ($y <= 1.5) {
				$res = self::logGamma1($y);
			} elseif ($y <= 4.0) {
				$res = self::logGamma2($y);
			} elseif ($y <= 12.0) {
				$res = self::logGamma3($y);
			} else {
				$res = self::logGamma4($y);
			}
		} else {
			// --------------------------
			//    Return for bad arguments
			// --------------------------
			$res = self::MAX_VALUE;
		}

		// ------------------------------
		//    Final adjustments and return
		// ------------------------------
		self::$logGammaCacheX = $x;
		self::$logGammaCacheResult = $res;

		return $res;
	}

	private static function logGamma1(float $y): float
	{
		// ---------------------
		//    EPS .LT. X .LE. 1.5
		// ---------------------
		if ($y < self::PNT68) {
			$corr = -log($y);
			$xm1 = $y;
		} else {
			$corr = 0.0;
			$xm1 = $y - 1.0;
		}

		$xden = 1.0;
		$xnum = 0.0;
		if ($y <= 0.5 || $y >= self::PNT68) {
			for ($i = 0; $i < 8; ++$i) {
				$xnum = $xnum * $xm1 + self::LG_P1[$i];
				$xden = $xden * $xm1 + self::LG_Q1[$i];
			}

			return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden));
		}

		$xm2 = $y - 1.0;
		for ($i = 0; $i < 8; ++$i) {
			$xnum = $xnum * $xm2 + self::LG_P2[$i];
			$xden = $xden * $xm2 + self::LG_Q2[$i];
		}

		return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
	}

	private static function logGamma2(float $y): float
	{
		// ---------------------
		//    1.5 .LT. X .LE. 4.0
		// ---------------------
		$xm2 = $y - 2.0;
		$xden = 1.0;
		$xnum = 0.0;
		for ($i = 0; $i < 8; ++$i) {
			$xnum = $xnum * $xm2 + self::LG_P2[$i];
			$xden = $xden * $xm2 + self::LG_Q2[$i];
		}

		return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
	}

	protected static function logGamma3(float $y): float
	{
		// ----------------------
		//    4.0 .LT. X .LE. 12.0
		// ----------------------
		$xm4 = $y - 4.0;
		$xden = -1.0;
		$xnum = 0.0;
		for ($i = 0; $i < 8; ++$i) {
			$xnum = $xnum * $xm4 + self::LG_P4[$i];
			$xden = $xden * $xm4 + self::LG_Q4[$i];
		}

		return self::LG_D4 + $xm4 * ($xnum / $xden);
	}

	protected static function logGamma4(float $y): float
	{
		// ---------------------------------
		//    Evaluate for argument .GE. 12.0
		// ---------------------------------
		$res = 0.0;
		if ($y <= self::LG_FRTBIG) {
			$res = self::LG_C[6];
			$ysq = $y * $y;
			for ($i = 0; $i < 6; ++$i) {
				$res = $res / $ysq + self::LG_C[$i];
			}
			$res /= $y;
			$corr = log($y);
			$res = $res + log(self::SQRT2PI) - 0.5 * $corr;
			$res += $y * ($corr - 1.0);
		}

		return $res;
	}
}

```