Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
// Apply Newton-Raphson step
$result = self::calculateDistribution($x, $alpha, $beta, true);
if (!is_float($result)) {
return ExcelError::NA();
}
$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 (!is_float($pdf)) {
return ExcelError::NA();
}
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
{
static $max = 32;
$summer = 0;
for ($n = 0; $n <= $max; ++$n) {
$divisor = $a;
for ($i = 1; $i <= $n; ++$i) {
$divisor *= ($a + $i);
}
$summer += ($x ** $n / $divisor);
}
return $x ** $a * exp(0 - $x) * $summer;
}
//
// Implementation of the Gamma function
//
public static function gammaValue(float $value): float
{
if ($value == 0.0) {
return 0;
}
static $p0 = 1.000000000190015;
static $p = [
1 => 76.18009172947146,
2 => -86.50532032941677,
3 => 24.01409824083091,
4 => -1.231739572450155,
5 => 1.208650973866179e-3,
6 => -5.395239384953e-6,
];
$y = $x = $value;
$tmp = $x + 5.5;
$tmp -= ($x + 0.5) * log($tmp);
$summer = $p0;
for ($j = 1; $j <= 6; ++$j) {
$summer += ($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
/** @var float */
private static $logGammaCacheResult = 0.0;
/** @var float */
private static $logGammaCacheX = 0.0;
/**
* logGamma function.
*
* Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
*
* The natural logarithm of the gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz
* Applied Mathematics Division
* Argonne National Laboratory
* Argonne, IL 60439
*
* References: *
* From the original documentation: *
** 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. *
*
* Error returns:
* 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.
*