Inconsistent treatment of s=0,∞ by gammaIncomplete(s, x)

[tedgin]

The regularized lower incomplete gamma function P(s,x) is defined as 𝛾(s,x)/𝛤(s) where x ∊ ℝ and s ∊ ℂ with re(s) > 0. 𝛾 is the lower incomplete gamma function, and 𝛤 is the gamma function.

According to its documentation, std.mathspecial.gammaIncomplete(s, x) implements P(s,x) where s and x have been restricted to ℝ⁺. However, it supports x = 0 and x = ∞, where gammaIncomplete(s, real.infinity) evaluates to lim_x→∞ P(s,x) = 1 for finite s. For the shape parameter s, its support for s = 0 and s = ∞ is inconsistent. When s is one of these values, for some values of x, gammaIncomplete(s, x) evaluates to NaN, and for others, it evaluates to a number.

P(s,x) = gammaIncomplete(s, x)

P(0⁺,0) = 0
P(0⁺,.1) = -nan
P(0⁺,1) = -nan
P(0⁺,10) = 1
P(0⁺,∞) = 1

P(∞,0) = 0
P(∞,.1) = 0
P(∞,1) = -nan
P(∞,10) = -nan
P(∞,∞) = -nan

Since s must be positive according to the definition of P, P(0,x) is undefined. However, as long as x > 0, lim_s→0⁺ P(s,x) = 1. When x = 0, lim_s→0⁺ P(s,0) = 0. Also, lim_s→∞ P(s,x) = 0, when x is non-negative and finite. lim_s,x→∞ P(s,x) does not exist, but lim_s→∞ lim_x→∞ P(s,x) = 1. To support using P as the gamma distribution CDF F(x; 𝛼,𝜆) ≝ P(𝛼,𝜆x) which requires that F(∞; 𝛼,𝜆) = 1, P(∞,∞) is often defined as lim_s→∞ lim_x→∞ P(s,x).

This issue proposes three changes to gammaIncomplete.

1. It is modified so that gammaIncomplete(+0., x) evaluates to 0 when x == 0 and 1 when x > 0.
2. It is modified so that gammaIncomplete(real.infinity, x) evaluates to 0, when x < real.infinity and 1 when x is real.infinity.
3. Its documentation is updated to reflect its existing support for x = 0.