Fix #10940 - Add support for degenerate cases to gammaIncomplete and gammaIncompleteCompl

std/internal/math/gammafunction.d

@@ -22,6 +22,9 @@ module std.internal.math.gammafunction;
 import std.internal.math.errorfunction;
 import std.math;
 import core.math : fabs, sin, sqrt;
+import std.algorithm.iteration : fold;
+import std.algorithm.searching : any;
+import std.range : only;
 
 pure:
 nothrow:
@@ -77,6 +80,34 @@ immutable real[8] logGammaDenominator = [
     -0x1.00f95ced9e5f54eep+9L, 1.0L
 ];
 
+
+/* Given a set of real values where at least one of them is NaN, it returns the
+ * NaN with the largest payload. This implements the NaN handling policy
+ * followed by D operators that accept multiple floating point arguments.
+ *
+ * It preserves the sign of the NaN. Also, when multiple NaNs have the largest
+ * payload, it returns the leftmost one. This means that
+ * largestNaNPayload(-NaN(1), NaN(1)) returns -NaN(1), and
+ * largestNaNPayload(NaN(1), -NaN(1)) returns NaN(1).
+ *
+ * When none of the provided values are NaN, the result is undefined.
+ */
+real largestNaNPayload(real first, real[] rest ...)
+{
+    return fold!((a, b) => cmp(abs(a), abs(b)) >= 0 ? a : b)(rest, first);
+}
+
+@safe unittest
+{
+    assert(largestNaNPayload(NaN(0)) is NaN(0));
+    assert(largestNaNPayload(NaN(1), NaN(0)) is NaN(1));
+    assert(largestNaNPayload(NaN(2), NaN(3), NaN(1)) is NaN(3));
+    assert(largestNaNPayload(-10.0L, -real.nan, 1.0L) is -real.nan);
+    assert(largestNaNPayload(-NaN(1), NaN(1)) is -NaN(1));
+    assert(largestNaNPayload(NaN(1), -NaN(1)) is NaN(1));
+}
+
+
 /*
  * Helper function: Gamma function computed by Stirling's formula.
  *
@@ -1012,14 +1043,7 @@ enum real BETA_BIGINV = 1.084202172485504434007e-19L;
  */
 real betaIncomplete(real aa, real bb, real xx )
 {
-    // If any parameters are NaN, return the NaN with the largest payload.
-    if (isNaN(aa) || isNaN(bb) || isNaN(xx))
-    {
-        // With cmp,
-        // -NaN(larger) < -NaN(smaller) < -inf < inf < NaN(smaller) < NaN(larger).
-        const largerParam = cmp(abs(aa), abs(bb)) >= 0 ? aa : bb;
-        return cmp(abs(xx), abs(largerParam)) >= 0 ? xx : largerParam;
-    }
+    if (only(aa, bb, xx).any!isNaN) return largestNaNPayload(xx, aa, bb);
 
     // domain errors
     if (signbit(aa) == 1 || signbit(bb) == 1) return real.nan;
@@ -1762,37 +1786,33 @@ real betaDistPowerSeries(real a, real b, real x )
 /***************************************
  *  Incomplete gamma integral and its complement
  *
- * These functions are defined by
- *
- *   gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
- *
- *  gammaIncompleteCompl(a,x)   =   1 - gammaIncomplete(a,x)
- * = ($(INTEGRATE x, &infin;) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
- *
- * In this implementation both arguments must be positive.
  * The integral is evaluated by either a power series or
  * continued fraction expansion, depending on the relative
  * values of a and x.
  */
 real gammaIncomplete(real a, real x )
 in
 {
-   assert(x >= 0);
-   assert(a > 0);
+    if (!any!isNaN(only(a, x)))
+    {
+        assert(x >= 0);
+        assert(signbit(a) == 0);
+    }
 }
 do
 {
-    /* left tail of incomplete gamma function:
-     *
-     *          inf.      k
-     *   a  -x   -       x
-     *  x  e     >   ----------
-     *           -     -
-     *          k=0   | (a+k+1)
-     *
-     */
-    if (x == 0)
-       return 0.0L;
+    // pass x first, so that if x and a are NaNs with the same payload but with
+    // opposite signs, return x.
+    if (any!isNaN(only(a, x))) return largestNaNPayload(x, a);
+
+    // domain violation
+    if (signbit(a) == 1 || x < 0.0L) return real.nan;
+
+    // edge cases
+    if (x == 0.0L) return 0.0L;
+    if (x is real.infinity) return 1.0L;
+    if (a is +0.0L) return 1.0L;
+    if (a is real.infinity) return 0.0L;
 
     if ( (x > 1.0L) && (x > a ) )
         return 1.0L - gammaIncompleteCompl(a,x);
@@ -1823,13 +1843,27 @@ do
 real gammaIncompleteCompl(real a, real x )
 in
 {
-   assert(x >= 0);
-   assert(a > 0);
+    if (!any!isNaN(only(a, x)))
+    {
+        assert(x >= 0);
+        assert(signbit(a) == 0);
+    }
 }
 do
 {
-    if (x == 0)
-        return 1.0L;
+    // pass x first, so that if x and a are NaNs with the same payload but with
+    // opposite signs, return x.
+    if (any!isNaN(only(a, x))) return largestNaNPayload(x, a);
+
+    // domain violation
+    if (signbit(a) == 1 || x < 0.0L) return real.nan;
+
+    // edge cases
+    if (x == 0.0L) return 1.0L;
+    if (x is real.infinity) return 0.0L;
+    if (a is +0.0L) return 0.0L;
+    if (a is real.infinity) return 1.0L;
+
     if ( (x < 1.0L) || (x < a) )
         return 1.0L - gammaIncomplete(a,x);
 
@@ -2039,6 +2073,18 @@ ihalve:
 
 @safe unittest
 {
+assert(gammaIncomplete(NaN(4), -NaN(4)) is -NaN(4));
+assert(!isNaN(gammaIncomplete(+0.0L, 1.0L)));
+assert(!isNaN(gammaIncomplete(1.0L, -0.0L)));
+assert(gammaIncomplete(+0.0L, 0.0L) == 0.0L);
+assert(gammaIncomplete(real.infinity, real.infinity) == 1.0L);
+
+assert(gammaIncompleteCompl(NaN(4), -NaN(4)) is -NaN(4));
+assert(!isNaN(gammaIncompleteCompl(+0.0L, 1.0L)));
+assert(!isNaN(gammaIncompleteCompl(1.0L, -0.0L)));
+assert(gammaIncompleteCompl(+0.0L, 0.0L) == 1.0L);
+assert(gammaIncompleteCompl(real.infinity, real.infinity) == 0.0L);
+
 //Values from Excel's GammaInv(1-p, x, 1)
 assert(fabs(gammaIncompleteComplInv(1, 0.5L) - 0.693147188044814L) < 0.00000005L);
 assert(fabs(gammaIncompleteComplInv(12, 0.99L) - 5.42818075054289L) < 0.00000005L);

std/mathspecial.d

@@ -29,6 +29,7 @@
  *      NAN = $(RED NAN)
  *      SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
  *      GAMMA = &#915;
+ *      LGAMMA =&#947;
  *      PSI = &Psi;
  *      THETA = &theta;
  *      INTEGRAL = &#8747;
@@ -47,9 +48,11 @@
  *      LT = &lt;
  *      LE = &le;
  *      GT = &gt;
+ *      GE = &ge;
  *      SQRT = &radic;
  *      HALF = &frac12;
  *      COMPLEX = &#8450;
+ *      NOBR = <nobr>$1</nobr>
  *
  * Copyright: Based on the CEPHES math library, which is
  *            Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
@@ -396,43 +399,112 @@ real betaIncompleteInverse(real a, real b, real y )
     return std.internal.math.gammafunction.betaIncompleteInv(a, b, y);
 }
 
-/** Incomplete gamma integral and its complement
+/** Regularized lower incomplete gamma function P(a,x)
  *
- * These functions are defined by
+ * Mathematically, P(a,x) = $(LGAMMA)(a,x)/$(GAMMA)(a), where $(LGAMMA)(a,x) is the lower incomplete
+ * gamma function, $(LGAMMA)(a,x) = $(INTEGRATE 0, x)$(POWER t, a-1)$(POWER e, -t)dt, a $(GT) 0, and
+ * x $(GE) 0.
  *
- *   gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
+ * Params:
+ *   a = the shape parameter, must be positive
+ *   x = the fraction of integration completion from below, must be non-negative
+ *
+ * Returns:
+ *   It returns P(a,x), an element of [0,1].
  *
- *  gammaIncompleteCompl(a,x)   =   1 - gammaIncomplete(a,x)
- * = ($(INTEGRATE x, $(INFIN)) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
+ * $(TABLE_SV
+ *   $(TR $(TH a)         $(TH x)              $(TH gammaIncomplete(a, x)) )
+ *   $(TR $(TD negative)  $(TD)                $(TD $(NAN))                )
+ *   $(TR $(TD)           $(TD $(LT) 0)        $(TD $(NAN))                )
+ *   $(TR $(TD positive)  $(TD 0)              $(TD 0)                     )
+ *   $(TR $(TD positive)  $(TD $(INFIN))       $(TD 1)                     )
+ *   $(TR $(TD +0)        $(TD $(GT) 0)        $(TD 1)                     )
+ *   $(TR $(TD $(INFIN))  $(TD (0, $(INFIN)))  $(TD 0)                     )
+ * )
  *
- * In this implementation both arguments must be positive.
- * The integral is evaluated by either a power series or
- * continued fraction expansion, depending on the relative
- * values of a and x.
+ * See_Also: $(LREF gamma) and $(LREF gammaIncompleteCompl)
  */
 real gammaIncomplete(real a, real x )
 in
 {
-   assert(x >= 0);
-   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(x))
+    {
+        assert(signbit(a) == 0, "a must be positive");
+        assert(x >= 0, "x must be non-negative");
+    }
 }
+out(p; isNaN(p) || (p >= 0 && p <= 1))
 do
 {
     return std.internal.math.gammafunction.gammaIncomplete(a, x);
 }
 
-/** ditto */
-real gammaIncompleteCompl(real a, real x )
+///
+@safe unittest
+{
+    assert(gammaIncomplete(1, 1) == 1 - 1/E);
+    assert(gammaIncomplete(1, 0) == 0);
+    assert(gammaIncomplete(1, real.infinity) == 1);
+    assert(gammaIncomplete(+0., 1) == 1);
+    assert(gammaIncomplete(real.infinity, 1) == 0);
+}
+
+/** Regularized upper incomplete gamma function Q(a,x)
+ *
+ * Mathematically, Q(a,x) = $(GAMMA)(a,x)/$(GAMMA)(a), where $(GAMMA)(a,x) is the upper incomplete
+ * gamma function, $(GAMMA)(a,x) = $(INTEGRATE x, $(INFIN))$(POWER t, a-1)$(POWER e, -t)dt,
+ * $(NOBR a $(GT) 0), and x $(GE) 0. Note that P(a,x) + Q(a,x) = 1 or Q(a,x) = 1 - P(a,x), so Q is
+ * the complement of P.
+ *
+ * Params:
+ *   a = the shape parameter, must be positive
+ *   x = the fraction of integration completion from above, must be non-negative
+ *
+ * Returns:
+ *   It returns Q(a,x), an element of [0,1].
+ *
+ * $(TABLE_SV
+ *   $(TR $(TH a)         $(TH x)              $(TH gammaIncompleteCompl(a, x)) )
+ *   $(TR $(TD negative)  $(TD)                $(TD $(NAN))                     )
+ *   $(TR $(TD)           $(TD $(LT) 0)        $(TD $(NAN))                     )
+ *   $(TR $(TD positive)  $(TD 0)              $(TD 1)                          )
+ *   $(TR $(TD positive)  $(TD $(INFIN))       $(TD 0)                          )
+ *   $(TR $(TD +0)        $(TD $(GT) 0)        $(TD 0)                          )
+ *   $(TR $(TD $(INFIN))  $(TD (0, $(INFIN)))  $(TD 1)                          )
+ * )
+ *
+ * See_Also: $(LREF gamma) and $(LREF gammaIncomplete)
+ */
+ real gammaIncompleteCompl(real a, real x )
 in
 {
-   assert(x >= 0);
-   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(x))
+    {
+        assert(signbit(a) == 0, "a must be positive");
+        assert(x >= 0, "x must be non-negative");
+    }
 }
+out(q; isNaN(q) || (q >= 0 && q <= 1))
 do
 {
     return std.internal.math.gammafunction.gammaIncompleteCompl(a, x);
 }
 
+///
+@safe unittest
+{
+    assert(isClose(gammaIncompleteCompl(2, 1), 2/E));
+    assert(gammaIncompleteCompl(1, 0) == 1);
+    assert(gammaIncompleteCompl(1, real.infinity) == 0);
+    assert(gammaIncompleteCompl(+0., 1) == 0);
+    assert(gammaIncompleteCompl(real.infinity, 1) == 1);
+    assert(isClose(gammaIncompleteCompl(1, 2), 1-gammaIncomplete(1, 2)));
+}
+
 /** Inverse of complemented incomplete gamma integral
  *
  * Given a and p, the function finds x such that