Fix #11103: keep small lower-tail targets in betaIncompleteInv (#11117)

Interval halving switched to the complement (a and b swapped,
y0 = 1 - yy0) once x passed 0.95. For a small yy0 that subtraction
rounds (to exactly 1 for 1e-20), so the target was lost. Switch only
when 1 - yy0 is exact, or to undo an earlier switch, and cap the
unswitched search at 1 - real.epsilon as the switched one already is.

* Fix #11103: switch to the complement unless 1 - yy0 loses yy0

Blocking the switch for every yy0 <= 0.5 made the search stay near 1
without the complement, which lost precision with 64-bit reals:
betaIncompleteInv(1e10, 1, 0.01) gave 0.99999991898 instead of
0.99999999954. Switch as before, and skip it only when 1 - yy0 rounds
away yy0. Drop the 1 - real.epsilon cap: the search now returns
nextDown(1), which is closer, since the CDF there is still below the
probability in the coverage tests.

* Fix #11103: solve for 1 - x when 1 - yy0 loses yy0

Skipping the complement switch kept the search in x near 1, which still
lost precision: betaIncompleteInv(1e11, 1, 1e-6) with 64-bit reals and
(1e12, 1, 0.007) with 80-bit reals were off from the 7th digit. When
1 - yy0 is not exact, bisect t = 1 - x in [1 - x1, 1 - x0] instead, with
betaIncompleteComplSmallX for the CDF, so neither yy0 nor the CDF is
subtracted from 1. As suggested by jmh530.

Signed-off-by: Sarthak <sarthakshreshtha345@gmail.com>

Author: Sarthak-Shreshtha01

std/internal/math/gammafunction.d

@@ -1439,6 +1439,26 @@ private real betaSmallShapeComplement(real a, real b, real x)
     return -expm1(a * (logBx + gammaRatioCorrection - logGammaOverA + seriesOverA) - shift);
 }
 
+// Solve betaIncomplete(aa, bb, 1 - t) = yy0 for t in [tlo, thi] by bisection
+// and return 1 - t. The CDF near 1 comes from betaIncompleteComplSmallX, so
+// neither yy0 nor the CDF is ever subtracted from 1.
+private real betaIncompleteInvTail(real aa, real bb, real yy0, real tlo, real thi)
+{
+    for (int i = 0; i < 200; i++)
+    {
+        const t = tlo + 0.5L * (thi - tlo);
+        if (t <= tlo || t >= thi || 1.0L - tlo == 1.0L - thi)
+            break;
+        if (betaIncompleteComplSmallX(bb, aa, t) > yy0)
+            tlo = t;
+        else
+            thi = t;
+    }
+    const x = 1.0L - (tlo + 0.5L * (thi - tlo));
+    // The quantile of a probability below 1 is below 1.
+    return x < 1.0L ? x : nextDown(1.0L);
+}
+
 /** Inverse of incomplete beta integral
  *
  * Given y, the function finds x such that
@@ -1562,6 +1582,14 @@ ihalve:
             else
                 di = (y0 - y)/(yh - yl);
             dir += 1;
+            if ( x0 > 0.95L && rflg == 0
+                && fabs((1.0L - (1.0L - yy0)) - yy0) > 8 * real.epsilon * yy0 )
+            {
+                // 1 - yy0 has lost part of yy0 (all of it for yy0 = 1e-20),
+                // so do not switch to the complement: solve for 1 - x instead.
+                x = betaIncompleteInvTail(aa, bb, yy0, 1.0L - x1, 1.0L - x0);
+                goto done;
+            }
             if ( x0 > 0.95L )
             {
                 if ( rflg == 1 )
@@ -1759,6 +1787,16 @@ done:
     assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842L), 0.229121208190918L) >= real.mant_dig - 3);
     assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= real.mant_dig - 3);
 
+    // https://github.com/dlang/phobos/issues/11103
+    // small lower-tail probabilities whose quantile is close to 1
+    assert(isClose(betaIncompleteInv(1000, 1, 1e-20L), pow(1e-20L, 1.0L / 1000), 1e-15L));
+    assert(isClose(betaIncompleteInv(1000, 1, 1e-10L), pow(1e-10L, 1.0L / 1000), 1e-15L));
+    assert(isClose(betaIncomplete(1000, 2, betaIncompleteInv(1000, 2, 1e-20L)), 1e-20L, 1e-10L));
+    // a moderate probability whose quantile is very close to 1 (also with 64-bit reals)
+    assert(isClose(betaIncompleteInv(1e10L, 1, 0.01L), pow(0.01L, 1.0L / 1e10L), 1e-14L));
+    assert(isClose(betaIncompleteInv(1e11L, 1, 1e-6L), pow(1e-6L, 1.0L / 1e11L), 1e-15L));
+    assert(isClose(betaIncompleteInv(1e12L, 1, 0.007L), pow(0.007L, 1.0L / 1e12L), 1e-15L));
+
     // Coverage tests. I don't have correct values for these tests, but
     // these values cover most of the code, so they are useful for
     // regression testing.
@@ -1767,8 +1805,11 @@ done:
     // significant improvement over the original CEPHES code.
     static if (real.mant_dig == 64) // 80-bit reals
     {
-        assert(betaIncompleteInv(0.01L, 8e-48L, 5.45464e-20L) == 1-real.epsilon);
-        assert(betaIncompleteInv(0.01L, 8e-48L, 9e-26L) == 1-real.epsilon);
+        // The CDF at nextDown(1) is still below these probabilities, so the
+        // quantile lies above it and nextDown(1) is the closest result below 1.
+        assert(betaIncomplete(0.01L, 8e-48L, nextDown(1.0L)) < 9e-26L);
+        assert(betaIncompleteInv(0.01L, 8e-48L, 5.45464e-20L) == nextDown(1.0L));
+        assert(betaIncompleteInv(0.01L, 8e-48L, 9e-26L) == nextDown(1.0L));
 
         // Beware: a one-bit change in pow() changes almost all digits in the result!
         // scipy says that this is 0.99999_99995_89020_6 (0x1.ffff_fffc_783f_2a7ap-1)