Fix BETAINV false convergence when the Beta CDF underflows

calculateInverse() treated a CDF value of exactly 0.0 as "the guess is the
root" and collapsed the bracket with $b = $a, ending the search. inverse()
already rejects probability <= 0, so a CDF of 0.0 inside the bisection can
only be a float64 underflow at a guess far below the root - never an exact
hit. It is just an ordinary "guess too low" and belongs in the existing
else branch.

The underflow is reached whenever a probe lands many standard deviations
from the mean, so the search collapses after two or three iterations and
returns whichever midpoint it was holding. BETAINV(0.5, 5000, 5000) gave
0.25, though Beta(a, a) is symmetric and its median is exactly 0.5;
BETAINV(0.5, 20000, 3) gave 0.5 against a true 0.99986. Onset is around
alpha = 1080 for beta = 1, where the closed form 0.5 ** (1 / alpha) is
available to check against.

That test also happened to stop the search for shapes where incompleteBeta
declines to evaluate at all and returns 0 for every x, so inverse() now
rejects alpha + beta above that documented limit up front: the CDF is
identically zero there, so no quantile exists to search for.

Tests cover a shape x probability grid against scipy reference values plus
three checks that need no external oracle: the symmetric median, the
Beta(alpha, 1) closed form, and the BETADIST round trip / mirror identity.

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

@@ -116,6 +116,11 @@ class Beta
         if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0.0)) {
             return ExcelError::NAN();
         }
+        if (($alpha + $beta) > self::LOG_GAMMA_X_MAX_VALUE) {
+            // incompleteBeta declines to evaluate here and returns 0 for every x,
+            // so there is no quantile to search for.
+            return ExcelError::NAN();
+        }
 
         return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax);
     }
@@ -130,7 +135,7 @@ class Beta
         while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
             $guess = ($a + $b) / 2;
             $result = self::distribution($guess, $alpha, $beta);
-            if (($result === $probability) || ($result === 0.0)) {
+            if ($result === $probability) {
                 $b = $a;
             } elseif ($result > $probability) {
                 $b = $guess;

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

@@ -40,4 +40,76 @@ class BetaInvTest extends AllSetupTeardown
             ],
         ];
     }
+
+    /**
+     * Beta(alpha, alpha) is symmetric about 0.5, so its median is exactly 0.5.
+     */
+    #[\PHPUnit\Framework\Attributes\DataProvider('providerSymmetricShape')]
+    public function testBetaInvMedianOfSymmetricDistribution(float $alpha): void
+    {
+        $result = Calculation::getInstance()->calculateFormula("=BETAINV(0.5, $alpha, $alpha)");
+        self::assertEqualsWithDelta(0.5, $result, 1.0e-12);
+    }
+
+    public static function providerSymmetricShape(): array
+    {
+        $cases = [];
+        foreach ([0.5, 2, 20, 200, 1000, 2000, 5000, 10000, 20000, 100000] as $alpha) {
+            $cases["alpha = beta = $alpha"] = [(float) $alpha];
+        }
+
+        return $cases;
+    }
+
+    /**
+     * The Beta(alpha, 1) CDF is x ** alpha, so its inverse is probability ** (1 / alpha).
+     */
+    #[\PHPUnit\Framework\Attributes\DataProvider('providerClosedForm')]
+    public function testBetaInvClosedFormForBetaOne(float $probability, float $alpha): void
+    {
+        $result = Calculation::getInstance()->calculateFormula("=BETAINV($probability, $alpha, 1)");
+        self::assertEqualsWithDelta($probability ** (1 / $alpha), $result, 1.0e-12);
+    }
+
+    public static function providerClosedForm(): array
+    {
+        $cases = [];
+        foreach ([0.001, 0.05, 0.25, 0.5, 0.9, 0.999] as $probability) {
+            foreach ([1, 10, 500, 1100, 5000, 20000, 100000] as $alpha) {
+                $cases["probability $probability, alpha $alpha"] = [$probability, (float) $alpha];
+            }
+        }
+
+        return $cases;
+    }
+
+    /**
+     * BETADIST must map the quantile BETAINV produced back onto the probability
+     * it was asked for, and BETAINV(p, a, b) must mirror 1 - BETAINV(1 - p, b, a).
+     */
+    #[\PHPUnit\Framework\Attributes\DataProvider('providerShapeGrid')]
+    public function testBetaInvAgreesWithBetaDist(float $probability, float $alpha, float $beta): void
+    {
+        $calculation = Calculation::getInstance();
+
+        $roundTrip = $calculation->calculateFormula("=BETADIST(BETAINV($probability, $alpha, $beta), $alpha, $beta)");
+        self::assertEqualsWithDelta($probability, $roundTrip, 1.0e-6, 'round trip through BETADIST');
+
+        $mirror = $calculation->calculateFormula(
+            "=BETAINV($probability, $alpha, $beta) + BETAINV(" . (1 - $probability) . ", $beta, $alpha)"
+        );
+        self::assertEqualsWithDelta(1.0, $mirror, 1.0e-12, 'mirror identity');
+    }
+
+    public static function providerShapeGrid(): array
+    {
+        $cases = [];
+        foreach ([0.01, 0.1, 0.5, 0.9, 0.99] as $probability) {
+            foreach ([[0.5, 0.5], [2, 5], [50, 50], [1000, 2], [5000, 3], [5000, 5000], [20000, 50], [200, 20000]] as [$alpha, $beta]) {
+                $cases["probability $probability, alpha $alpha, beta $beta"] = [$probability, (float) $alpha, (float) $beta];
+            }
+        }
+
+        return $cases;
+    }
 }

tests/data/Calculation/Statistical/BETAINV.php

@@ -51,6 +51,60 @@ return [
         '#VALUE!',
         0.2, 4, 5, 0, 'NAN',
     ],
+    // Large shape parameters: the CDF underflows to 0 at the low probes of the
+    // bisection long before the root is reached. Onset is alpha ~ 1080 for beta = 1.
+    'symmetric, median' => [
+        0.5,
+        0.5, 5000, 5000, 0, 1,
+    ],
+    'symmetric, upper quartile' => [
+        0.503372494704,
+        0.75, 5000, 5000, 0, 1,
+    ],
+    'symmetric, lower decile' => [
+        0.493592345088,
+        0.1, 5000, 5000, 0, 1,
+    ],
+    'symmetric, median, alpha = beta = 20000' => [
+        0.5,
+        0.5, 20000, 20000, 0, 1,
+    ],
+    'symmetric, scaled to [1, 3]' => [
+        2.006744989408,
+        0.75, 5000, 5000, 1, 3,
+    ],
+    'beta = 1, closed form 0.5 ** (1 / alpha)' => [
+        0.999861380173,
+        0.5, 5000, 1, 0, 1,
+    ],
+    'beta = 1, just above the onset' => [
+        0.999370064692,
+        0.5, 1100, 1, 0, 1,
+    ],
+    'right-skewed, alpha >> beta' => [
+        0.999866312606,
+        0.5, 20000, 3, 0, 1,
+    ],
+    'left-skewed, beta >> alpha' => [
+        0.000133687394,
+        0.5, 3, 20000, 0, 1,
+    ],
+    'moderately large alpha' => [
+        0.973403556565,
+        0.25, 2000, 50, 0, 1,
+    ],
+    'alpha so large that the quantile rounds to 1' => [
+        1.0,
+        0.5, 1.0e20, 3, 0, 1,
+    ],
+    'huge symmetric shapes still give the median' => [
+        0.5,
+        0.5, 1.0e20, 1.0e20, 0, 1,
+    ],
+    'alpha + beta beyond the range incompleteBeta will evaluate' => [
+        '#NUM!',
+        0.5, 2.0e305, 1.0e305, 0, 1,
+    ],
     'alpha < 0' => [
         '#NUM!',
         0.2, -4, 5, 0, 1,