Expand GAMMAINV bracket so tail quantiles are not clamped

GammaBase::calculateInverse fixed its upper bound at alpha*beta*5, so any
GAMMA.INV/GAMMAINV quantile larger than that was clamped to it: e.g.
GAMMAINV(0.9999, 1, 1) returned 5 rather than ~9.2103, breaking the
round-trip GAMMADIST(GAMMAINV(p)) == p.

Grow the upper bound geometrically until it brackets the root. If the CDF
stops increasing first (the series approximation is past its usable range)
keep the original bound instead of expanding into it, which also stops a
probability the series cannot reach from running the bound away.

src/PhpSpreadsheet/Calculation/Statistical/Distributions/GammaBase.php

@@ -32,6 +32,24 @@ abstract class GammaBase
         $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;

tests/PhpSpreadsheetTests/Calculation/Functions/Statistical/GammaInvTest.php

@@ -5,9 +5,48 @@ declare(strict_types=1);
 namespace PhpOffice\PhpSpreadsheetTests\Calculation\Functions\Statistical;
 
 use PhpOffice\PhpSpreadsheet\Calculation\Calculation;
+use PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions\Gamma;
 
 class GammaInvTest extends AllSetupTeardown
 {
+    /**
+     * Extreme upper-tail quantiles whose true root exceeds the old fixed
+     * alpha*beta*5 bracket ceiling, which used to clamp the result to it.
+     * Expected values from mpmath (findroot on the regularized gammainc).
+     */
+    #[\PHPUnit\Framework\Attributes\DataProvider('providerGammaInvExtremeTail')]
+    public function testGammaInvExtremeTail(float $expected, float $probability, float $alpha, float $beta): void
+    {
+        $x = Gamma::inverse($probability, $alpha, $beta);
+        self::assertIsFloat($x);
+        // Bracket ceiling was exceeded, so the fix must have expanded past it.
+        self::assertGreaterThan($alpha * $beta * 5.0, $x);
+        // Round-trip invariant: the quantile maps back to the input probability.
+        $roundTrip = Gamma::distribution($x, $alpha, $beta, true);
+        self::assertEqualsWithDelta($probability, $roundTrip, 1.0e-8);
+        // And it matches the reference quantile.
+        self::assertEqualsWithDelta($expected, $x, 1.0e-3);
+    }
+
+    public static function providerGammaInvExtremeTail(): array
+    {
+        return [
+            'p=0.9999 alpha=1 beta=1' => [9.210340371976, 0.9999, 1.0, 1.0],
+            'p=0.99995 alpha=1 beta=1' => [9.903487552536, 0.99995, 1.0, 1.0],
+            'p=0.9999 alpha=0.5 beta=2' => [15.136705226623, 0.9999, 0.5, 2.0],
+            'p=0.9999 alpha=1 beta=2' => [18.420680743952, 0.9999, 1.0, 2.0],
+        ];
+    }
+
+    public function testGammaInvUnreachableTailStaysBounded(): void
+    {
+        // A probability the forward series cannot reach must not send the
+        // bracket expansion running away to a huge nonsensical quantile.
+        $x = Gamma::inverse(0.9999999, 1.0, 1.0);
+        self::assertIsFloat($x);
+        self::assertLessThan(1000.0, $x);
+    }
+
     #[\PHPUnit\Framework\Attributes\DataProvider('providerGAMMAINV')]
     public function testGAMMAINV(mixed $expectedResult, mixed ...$args): void
     {