upper incomplete gamma function Γ(a,z) (#271)

* clarify gamma_inc docs and make ind=0 the default, rm incorrect method that converted all Integer (even BigInt) to Float64

* add gamma(a,x)

* more promotion fixes

* gamma docs and cross-refs

* more promotion tweaks

* use inc=0 default parameter in gamma_inc tests

* more fixes and tests

* more tests

* doc typo

* rm unused method

* add more tests and infinity handling

* whoops, move tests out of loop

Author: stevengj

docs/src/functions_list.md

@@ -55,7 +55,8 @@ SpecialFunctions.ellipk
 SpecialFunctions.ellipe
 SpecialFunctions.eta
 SpecialFunctions.zeta
-SpecialFunctions.gamma
+SpecialFunctions.gamma(::Number)
+SpecialFunctions.gamma(::Number,::Number)
 SpecialFunctions.gamma_inc
 SpecialFunctions.gamma_inc_inv
 SpecialFunctions.beta_inc

docs/src/functions_overview.md

@@ -5,11 +5,12 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 ## [Gamma Function](https://dlmf.nist.gov/5)
 | Function | Description |
 |:-------- |:----------- |
-| [`gamma(z)`](@ref SpecialFunctions.gamma) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` |
+| [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` |
 | [`digamma(x)`](@ref SpecialFunctions.digamma)                 | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `lgamma` at `x`)                                                     |
 | [`invdigamma(x)`](@ref SpecialFunctions.invdigamma)   | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) |
 | [`trigamma(x)`](@ref SpecialFunctions.trigamma)     | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) |
 | [`polygamma(m,x)`](@ref SpecialFunctions.polygamma)  | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `lgamma` function at `x`) |
+| [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number))  | [upper incomplete gamma function ``\Gamma(a,z)``](https://en.wikipedia.org/wiki/Incomplete_gamma_function) |
 | [`gamma_inc(a,x,IND)`](@ref SpecialFunctions.gamma_inc)  | [incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates P(a,x) and Q(a,x)for accuracy specified by IND and returns tuple (p,q)) |
 | [`beta_inc(a,b,x,y)`](@ref SpecialFunctions.beta_inc)  | [incomplete beta function ratio Ix(a,b) and Iy(a,b)](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q)) |
 | [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv)  | [inverse of incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates x given P(a,x)=p and Q(a,x)=q  |

src/gamma.jl

@@ -573,7 +573,6 @@ export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, lo
 gamma(x::Float64)       = nan_dom_err(ccall((:tgamma, libm), Float64, (Float64,), x), x)
 gamma(x::Float32)       = nan_dom_err(ccall((:tgammaf, libm), Float32, (Float32,), x), x)
 gamma(x::Float16)       = Float16(gamma(Float32(x)))
-gamma(x::Real)          = gamma(float(x))
 gamma(x::AbstractFloat) = throw(MethodError(gamma, x))
 
 function gamma(x::BigFloat)
@@ -606,15 +605,17 @@ External links:
 [DLMF](https://dlmf.nist.gov/5.2.1),
 [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function).
 
-See also: [`loggamma(z)`](@ref SpecialFunctions.loggamma).
+See also: [`loggamma(z)`](@ref SpecialFunctions.loggamma) for ``\log \Gamma(z)`` and
+[`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number)) for
+the upper incomplete gamma function ``\\Gamma(a,z)``.
 
 # Implementation by
 - `Float`: C standard math library
     [libm](https://en.wikipedia.org/wiki/C_mathematical_functions#libm).
 - `Complex`: by `exp(loggamma(z))`.
 - `BigFloat`: C library for multiple-precision floating-point [MPFR](https://www.mpfr.org/)
 """
-gamma
+gamma(x::Number) = gamma(float(x))
 
 function logabsgamma(x::Float64)
     signp = Ref{Int32}()

src/gamma_inc.jl

@@ -818,7 +818,7 @@ end
 # doi>10.1145/22721.23109
 
 """
-    gamma_inc(a,x,IND)
+    gamma_inc(a,x,IND=0)
 
 Returns a tuple ``(p, q)`` where ``p + q = 1``, and
 ``p=P(a,x)`` is the Incomplete gamma function ratio given by:
@@ -829,14 +829,19 @@ and ``q=Q(a,x)`` is the Incomplete gamma function ratio given by:
 ```math
 Q(x,a)=\\frac{1}{\\Gamma (a)} \\int_{x}^{\\infty} e^{-t}t^{a-1} dt.
 ```
+In terms of these, the lower incomplete gamma function is
+``\\gamma(a,x) = P(a,x) \\Gamma(a)`` and the upper incomplete
+gamma function is ``\\Gamma(a,x) = Q(a,x) \\Gamma(a)``.
 
 `IND ∈ [0,1,2]` sets accuracy: `IND=0` means 14 significant digits accuracy, `IND=1` means 6 significant digit, and `IND=2` means only 3 digit accuracy.
 
 External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)
 
 See also [`gamma(z)`](@ref SpecialFunctions.gamma), [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv)
 """
-function gamma_inc(a::Float64,x::Float64,ind::Integer)
+gamma_inc(a::Real,x::Real,ind::Integer=0) = _gamma_inc(promote(float(a),float(x))...,ind)
+
+function _gamma_inc(a::Float64,x::Float64,ind::Integer)
     iop = ind + 1
     acc = acc0[iop]
     if a<0.0 || x<0.0
@@ -933,17 +938,14 @@ function gamma_inc(a::Float64,x::Float64,ind::Integer)
 
 end
 
-function gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall
+function _gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall
     z = BigFloat()
     ccall((:mpfr_gamma_inc, :libmpfr), Int32 , (Ref{BigFloat} , Ref{BigFloat} , Ref{BigFloat} , Int32) , z , a , x , ROUNDING_MODE[])
     q = z/gamma(a)
     return (1.0 - q, q)
 end
-gamma_inc(a::Float32,x::Float32,ind::Integer) = ( Float32(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float32(gamma_inc(Float64(a),Float64(x),ind)[2]) )
-gamma_inc(a::Float16,x::Float16,ind::Integer) = ( Float16(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float16(gamma_inc(Float64(a),Float64(x),ind)[2]) )
-gamma_inc(a::Real,x::Real,ind::Integer) = (gamma_inc(float(a),float(x),ind))
-gamma_inc(a::Integer,x::Integer,ind::Integer) = gamma_inc(Float64(a),Float64(x),ind)
-gamma_inc(a::AbstractFloat,x::AbstractFloat,ind::Integer) = throw(MethodError(gamma_inc,(a,x,ind,"")))
+_gamma_inc(a::Float32,x::Float32,ind::Integer) = Float32.(_gamma_inc(Float64(a),Float64(x),ind))
+_gamma_inc(a::Float16,x::Float16,ind::Integer) = Float16.(_gamma_inc(Float64(a),Float64(x),ind))
 
 #EFFICIENT AND ACCURATE ALGORITHMS FOR THECOMPUTATION AND INVERSION OF THE INCOMPLETEGAMMA FUNCTION RATIOS by Amparo Gil, Javier Segura, Nico M. Temme
 #SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981
@@ -958,7 +960,9 @@ External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wiki
 
 See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc).
 """
-function gamma_inc_inv(a::Float64, p::Float64, q::Float64)
+gamma_inc_inv(a::Real, p::Real, q::Real) = _gamma_inc_inv(promote(float(a), float(p), float(q))...)
+
+function _gamma_inc_inv(a::Float64, p::Float64, q::Float64)
     if p < 0.5
         pcase = true
         porq = p
@@ -1028,8 +1032,70 @@ function gamma_inc_inv(a::Float64, p::Float64, q::Float64)
     return x
 end
 
-gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
-gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
-gamma_inc_inv(a::Real, p::Real, q::Real) = ( gamma_inc_inv(float(a), float(p), float(q)) )
-gamma_inc_inv(a::Integer, p::Integer, q::Integer) = ( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
-gamma_inc_inv(a::AbstractFloat, p::AbstractFloat, q::AbstractFloat) = throw(MethodError(gamma_inc_inv,(a,p,q,"")))
+_gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
+_gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16(_gamma_inc_inv(Float64(a), Float64(p), Float64(q)))
+
+# like promote(x,y), but don't complexify real values
+promotereal(x::Real,y::Real) = promote(x,y)
+promotereal(x::Real,y::Number) = let (u,v) = promote(x,y); real(u),v; end
+promotereal(x::Number,y::Real) = let (u,v) = promote(x,y); u,real(v); end
+promotereal(x,y) = promote(x,y)
+
+"""
+    gamma(a,x)
+
+Returns the upper incomplete gamma function
+```math
+\\Gamma(a,x) = \\int_x^\\infty t^{a-1} e^{-t} dt \\, ,
+```
+supporting arbitrary real or complex `a` and `x`.
+
+(The ordinary gamma function [`gamma(x)`](@ref) corresponds to ``\\Gamma(a) = \\Gamma(a,0)``.
+See also the [`gamma_inc`](@ref) function to compute both the upper and lower
+(``\\gamma(a,x)``) incomplete gamma functions scaled by ``\\Gamma(a)``.
+
+External links: [DLMF](https://dlmf.nist.gov/8.2.2), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)
+"""
+gamma(a::Number,x::Number) = _gamma(promotereal(float(a),float(x))...)
+gamma(a::Integer,x::Number) = _gamma(a, float(x))
+
+function _gamma(a::Number,x::Number)
+    if a isa Real && x isa Real && !isfinite(a*x)
+        if isinf(x) && isfinite(a)
+            if x > 0 # == +Inf
+                return one(a)*zero(x)
+            elseif a > 0 && isinteger(a) # x == -Inf
+                return one(a)*x # -Inf
+            end
+        elseif isinf(a) && isfinite(x)
+            if a > 0
+                if x ≥ 0
+                    return a*one(x) # +Inf
+                else
+                    throw(DomainError((a,x), "gamma will only return a complex result if called with a complex argument"))
+                end
+            elseif a < 0
+                return zero(a)*one(x)
+            end
+        end
+    end
+    return iszero(x) ? one(x)*gamma(a) : x^a * expint(1-a, x)
+end
+
+_gamma(a::Integer, x::BigFloat) = _gamma_big(a, x)
+_gamma(a::BigInt, x::Real) = _gamma_big(a, x)
+_gamma(a::BigInt, x::BigFloat) = _gamma_big(a, x)
+_gamma(a::BigFloat,x::BigFloat) = _gamma_big(a, x)
+function _gamma_big(a::Real,x::Real)
+    if x < 0
+        # MPFR returns NaN in this case
+        if isinteger(a) && a > 0
+            return invoke(_gamma, Tuple{Number,Number}, a, BigFloat(x))
+        else
+            throw(DomainError((a,x), "gamma will only return a complex result if called with a complex argument"))
+        end
+    end
+    z = BigFloat()
+    ccall((:mpfr_gamma_inc, :libmpfr), Int32 , (Ref{BigFloat} , Ref{BigFloat} , Ref{BigFloat} , Int32) , z , a , x , ROUNDING_MODE[])
+    return z
+end

test/gamma_inc.jl

@@ -1,41 +1,41 @@
 @testset "incomplete gamma ratios" begin
 #Computed using Wolframalpha gamma(a,x)/gamma(a) ~ gamma_q(a,x,0) function.
-    @test gamma_inc(10,10,0)[2] ≈ 0.45792971447185221
-    @test gamma_inc(1,1,0)[2] ≈ 0.3678794411714423216
-    @test gamma_inc(0.5,0.5,0)[2] ≈ 0.31731050786291410
-    @test gamma_inc(BigFloat(30.5),BigFloat(30.5),0)[2] ≈ parse(BigFloat,"0.47591691193354987004") rtol=eps()
-    @test gamma_inc(5.5,0.5,0)[2] ≈ 0.9999496100513121669
-    @test gamma_inc(0.5,7.4,0)[2] ≈ 0.0001195355018130302
-    @test gamma_inc(0.5,0.22,0)[2] ≈ 0.507122455359825146
-    @test gamma_inc(0.5,0.8,0)[2] ≈ 0.20590321073206830887
-    @test gamma_inc(11.5,0.5,0)[2] ≈ 0.999999999998406112
-    @test gamma_inc(0.19,0.99,0)[2] ≈ 0.050147247342905857
-    @test gamma_inc(0.9999,0.9999,0)[2] ≈ 0.3678730556923103
-    @test gamma_inc(24,23.9999999999,0)[2] ≈ 0.472849720555859138
-    @test gamma_inc(0.5,0.55,0)[2] ≈ 0.29426610430496289
-    @test gamma_inc(Float32(0.5),Float32(0.55),0)[2] ≈ Float32(gamma_inc(0.5,0.55,0)[2])
-    @test gamma_inc(Float16(0.5),Float16(0.55),0)[2] ≈ Float16(gamma_inc(0.5,0.55,0)[2])
-    @test gamma_inc(30,29.99999,0)[2] ≈ 0.475717712451705704
-    @test gamma_inc(30,29.9,0)[2] ≈ 0.482992166284958565
-    @test gamma_inc(10,0.0001,0)[2] ≈ 1.0000
-    @test gamma_inc(0.0001,0.0001,0)[2] ≈ 0.000862958131006599
-    @test gamma_inc(0.0001,10.5,0)[1] ≈ 0.999999999758896146
-    @test gamma_inc(1,1,0)[1] ≈ 0.63212055882855768
-    @test gamma_inc(13,15.1,0)[2] ≈ 0.25940814264863701
-    @test gamma_inc(0.6,1.3,0)[2] ≈ 0.136458554006505355
-    @test gamma_inc((100),(80),0)[2] ≈ 0.9828916869648668
+    @test gamma_inc(10,10)[2] ≈ 0.45792971447185221
+    @test gamma_inc(1,1)[2] ≈ 0.3678794411714423216
+    @test gamma_inc(0.5,0.5)[2] ≈ 0.31731050786291410
+    @test gamma_inc(BigFloat(30.5),BigFloat(30.5))[2] ≈ parse(BigFloat,"0.47591691193354987004") rtol=eps()
+    @test gamma_inc(5.5,0.5)[2] ≈ 0.9999496100513121669
+    @test gamma_inc(0.5,7.4)[2] ≈ 0.0001195355018130302
+    @test gamma_inc(0.5,0.22)[2] ≈ 0.507122455359825146
+    @test gamma_inc(0.5,0.8)[2] ≈ 0.20590321073206830887
+    @test gamma_inc(11.5,0.5)[2] ≈ 0.999999999998406112
+    @test gamma_inc(0.19,0.99)[2] ≈ 0.050147247342905857
+    @test gamma_inc(0.9999,0.9999)[2] ≈ 0.3678730556923103
+    @test gamma_inc(24,23.9999999999)[2] ≈ 0.472849720555859138
+    @test gamma_inc(0.5,0.55)[2] ≈ 0.29426610430496289
+    @test gamma_inc(Float32(0.5),Float32(0.55))[2] ≈ Float32(gamma_inc(0.5,0.55)[2])
+    @test gamma_inc(Float16(0.5),Float16(0.55))[2] ≈ Float16(gamma_inc(0.5,0.55)[2])
+    @test gamma_inc(30,29.99999)[2] ≈ 0.475717712451705704
+    @test gamma_inc(30,29.9)[2] ≈ 0.482992166284958565
+    @test gamma_inc(10,0.0001)[2] ≈ 1.0000
+    @test gamma_inc(0.0001,0.0001)[2] ≈ 0.000862958131006599
+    @test gamma_inc(0.0001,10.5)[1] ≈ 0.999999999758896146
+    @test gamma_inc(1,1)[1] ≈ 0.63212055882855768
+    @test gamma_inc(13,15.1)[2] ≈ 0.25940814264863701
+    @test gamma_inc(0.6,1.3)[2] ≈ 0.136458554006505355
+    @test gamma_inc((100),(80))[2] ≈ 0.9828916869648668
     @test gamma_inc((100),(80),1)[2] ≈ 0.9828916869
     @test Float16(gamma_inc((100),(80),2)[2]) ≈ Float16(.983)
-    @test gamma_inc(13.5,15.1,0)[2] ≈ 0.305242642543419087
-    @test gamma_inc(11,9,0)[1] ≈ 0.2940116796594881834
-    @test gamma_inc(8,32,0)[1] ≈ 0.99999989060651042057
-    @test gamma_inc(15,16,0)[2] ≈ 0.3675273597655649298
-    @test gamma_inc(15.5,16,0)[2] ≈ 0.4167440299455427811
-    @test gamma_inc(0.9,0.8,0)[1] ≈ 0.59832030278768172
-    @test gamma_inc(1.7,2.5,0)[1] ≈ 0.78446115627678957
-    @test gamma_inc(11.1,0.001,0)[2] ≈ 1.0000
-    @test gamma_inc(1e7, (1e7)+1, 0)[1] ≈ 0.5001682088254367
-    @test gamma_inc(1e7, (1e7)+1, 0)[2] ≈ 0.4998317911745633
+    @test gamma_inc(13.5,15.1)[2] ≈ 0.305242642543419087
+    @test gamma_inc(11,9)[1] ≈ 0.2940116796594881834
+    @test gamma_inc(8,32)[1] ≈ 0.99999989060651042057
+    @test gamma_inc(15,16)[2] ≈ 0.3675273597655649298
+    @test gamma_inc(15.5,16)[2] ≈ 0.4167440299455427811
+    @test gamma_inc(0.9,0.8)[1] ≈ 0.59832030278768172
+    @test gamma_inc(1.7,2.5)[1] ≈ 0.78446115627678957
+    @test gamma_inc(11.1,0.001)[2] ≈ 1.0000
+    @test gamma_inc(1e7, (1e7)+1)[1] ≈ 0.5001682088254367
+    @test gamma_inc(1e7, (1e7)+1)[2] ≈ 0.4998317911745633
     @test_throws DomainError gamma_inc(-1,2,2)
     @test_throws DomainError gamma_inc(0,0,1)
 end
@@ -63,3 +63,40 @@ end
         @test SpecialFunctions.stirling(x) ≈ log(gamma(x)) - (x-.5)*log(x)+x- log(2*pi)/2.0
     end
 end
+
+double(x::Real) = Float64(x)
+double(x::Integer) = Int(x)
+double(x::Complex) = ComplexF64(x)
+
+@testset "upper incomplete gamma function" begin
+    setprecision(BigFloat, 256) do
+        for a in Any[0:0.4:3; 1:3], x in 0:0.2:2
+            @test gamma(a,x) ≈ gamma(big(a),big(x))
+        end
+        @test_throws DomainError gamma(-2.2, -1.3)
+        for (a,x, exact) in (
+            (2, big"+1.3", big"0.6268231239782289871813662853957039889809398944861850589869804057956189274569818"),
+            (2, big"-1.3", big"-1.100789000285773266137246974803445875429455665367265440176867948378982424804222"),
+            (big"2.2", big"1.3", big"0.7550934853735524106456916078787596171599416239996064979513848803118671763945577"),
+            (big"-2.2", big"1.3", big"0.03938564959393195337328006806473296774233286427565465045140458665248322893076784"),
+            (big"-2.2", big"-1.3"+0im, Complex(big"0.1688350167695890519002747177035271528716667324397453142169971691104641885370081", big"-1.724677939954414412829215929952389349929001266936564849801773346878175589599461")),
+            (2, big"1.3"+2im, Complex(big"0.2347744561498965806069446787000456827042253434143074251879541754828136789074293", big"-0.7967951407674886396960561446265346226800340382780520672369163777061700801594228")),
+            (big"2.2", big"1.3"+2im, Complex(big"0.4226969694490623767886715572749155659150206342283251276656790859161308520476838", big"-0.9507541450333763666696561254491461255244463810462729822382144498284499948680099")),
+            (big"2.2"-2im, big"1.3"+2im, Complex(big"-3.858698824275578861673404086439928163791926530897603434146764316225759116935376", big"0.2105970967650200831829566202893410569725743566037534990454186505024476246014339")),
+            (30, big"-1e2", big"-2.080142496354742431762923569133534773079637504483999408076147454085194055640495e+101"),
+            )
+            if !(exact isa Complex) # we don't yet have a Complex{BigFloat} gamma function
+                @test gamma(a, x) ≈ exact atol=eps(abs(exact))*1000
+            end
+            @test gamma(double(a), double(x)) ≈ double(exact) rtol=1e-13
+        end
+    end
+    @test gamma(30, big(-1000)) ≈ big"-1.914496653187781420578078670260160511219745144309147121829045802464973219500799e+521" rtol=1e-70
+    @test gamma(30, -1000) == -Inf
+    @test gamma(Inf, 51.2) == Inf
+    @test gamma(-Inf, -51.2) == 0.0
+    @test_throws DomainError gamma(Inf, -51.2)
+    @test gamma(2.3, Inf) == 0.0
+    @test gamma(2, -Inf) == -Inf
+    @test_throws DomainError gamma(2.2, -Inf)
+end