10962 improved handling of gammaIncompleteCompleInverse edge cases (dlang#10963)

* Made NaN handling conform to that of builtin binary operators, i.e.,
  when both input values are NaN, the one with the larger payload is
  returned.
* Extended domain coverage to handle a = +0
* Short circuit cases where Q isn't invertible and return NaN
* Explicitly handle p = 1 edge case.

Author: tedgin

std/internal/math/gammafunction.d

@@ -1945,12 +1945,29 @@ do
 real gammaIncompleteComplInv(real a, real p)
 in
 {
-  assert(p >= 0 && p <= 1);
-  assert(a>0);
+    if (!any!isNaN(only(a, p)))
+    {
+        assert(signbit(a) == 0);
+        assert(p >= 0.0L && p <= 1.0L);
+    }
 }
 do
 {
-    if (p == 0) return real.infinity;
+    // pass p first, so that if p and a are NaNs with the same payload but with
+    // opposite signs, return p.
+    if (any!isNaN(only(a, p))) return largestNaNPayload(p, a);
+
+    // domain violations
+    if (signbit(a) == 1) return real.nan;
+    if (p < 0.0L || p > 1.0L) return real.nan;
+
+    // places where not invertible
+    if (a is +0.0L && p < 1.0L) return real.nan;
+    if (a is real.infinity && p > 0.0L) return real.nan;
+
+    // edge cases for p
+    if (p == 0.0L) return real.infinity;
+    if (p == 1.0L) return 0.0L;
 
     real y0 = p;
     const real MAXLOGL =  1.1356523406294143949492E4L;
@@ -2107,6 +2124,13 @@ static if (real.mant_dig >= 64) // incl. 80-bit reals
     assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109L) < 0.000000000005L);
 else
     assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109L) < 0.00000005L);
+
+assert(gammaIncompleteComplInv(NaN(0x5UL), -NaN(0x5UL)) is -NaN(0x5UL));
+assert(!isNaN(gammaIncompleteComplInv(+0.0L, 1.0L)));
+assert(isNaN(gammaIncompleteComplInv(+0.0L, nextDown(1.0L))));
+assert(!isNaN(gammaIncompleteComplInv(real.infinity, -0.0L)));
+assert(isNaN(gammaIncompleteComplInv(real.infinity, nextUp(+0.0L))));
+assert(gammaIncompleteComplInv(2.0L, 1.0L) == 0.0L);
 }
 
 

std/mathspecial.d

@@ -505,23 +505,59 @@ do
     assert(isClose(gammaIncompleteCompl(1, 2), 1-gammaIncomplete(1, 2)));
 }
 
-/** Inverse of complemented incomplete gamma integral
+/** Inverse regularized upper incomplete gamma function Q$(SUP -1)(a,p) with respect to p
  *
- * Given a and p, the function finds x such that
+ * Given a and p, the function finds x such that p = Q(a,x).
  *
- *  gammaIncompleteCompl( a, x ) = p.
+ * Params:
+ *   a = the shape parameter, must be positive
+ *   p = Q(a,x), must be in the interval [0,1]
+ *
+ * Returns:
+ *   It returns x, a value $(GE) 0
+ *
+ * $(TABLE_SV
+ *   $(TR $(TH a)               $(TH p)        $(TH gammaIncompleteComplInverse(a, p)) )
+ *   $(TR $(TD negative)        $(TD)          $(TD $(NAN))                            )
+ *   $(TR $(TD)                 $(TD $(LT) 0)  $(TD $(NAN))                            )
+ *   $(TR $(TD)                 $(TD $(GT) 1)  $(TD $(NAN))                            )
+ *   $(TR $(TD +0)              $(TD $(LT) 1)  $(TD $(NAN))                            )
+ *   $(TR $(TD $(INFIN))        $(TD $(GT) 0)  $(TD $(NAN))                            )
+ *   $(TR $(TD $(GT) 0)         $(TD 0)        $(TD $(INFIN))                          )
+ *   $(TR $(TD $(LT) $(INFIN))  $(TD 1)        $(TD 0)                                 )
+ * )
+ *
+ * See_Also: $(LREF gammaIncompleteCompl)
  */
 real gammaIncompleteComplInverse(real a, real p)
 in
 {
-  assert(p >= 0 && p <= 1);
-  assert(a > 0);
+    // allow NaN input to pass through so that it can be addressed by the
+    // internal NaN payload propagation logic
+    if (!isNaN(a) && !isNaN(p))
+    {
+        assert(signbit(a) == 0, "a must be positive");
+        assert(p >= 0.0L && p <= 1.0L, "p must be in the interval [0,1]");
+    }
 }
+out(x; isNaN(x) || x >= 0.0L)
 do
 {
     return std.internal.math.gammafunction.gammaIncompleteComplInv(a, p);
 }
 
+///
+@safe unittest
+{
+    const a = 2, p = 0.5L;
+    assert(isClose(gammaIncompleteComplInverse(a, gammaIncompleteCompl(a, p)), p));
+
+    assert(gammaIncompleteComplInverse(1, 1/E) == 1);
+    assert(isNaN(gammaIncompleteComplInverse(+0.0L, 0.1)));
+    assert(isNaN(gammaIncompleteComplInverse(real.infinity, 0.2)));
+    assert(gammaIncompleteComplInverse(3, 0) is real.infinity);
+    assert(gammaIncompleteComplInverse(4, 1) == 0);
+}
 
 /* ***********************************************
  *     ERROR FUNCTIONS & NORMAL DISTRIBUTION     *