Increase test coverage, fix NaN

.github/workflows/CI.yml

@@ -36,7 +36,7 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
           show-versioninfo: true
-      - uses: julia-actions/cache@v2
+      - uses: julia-actions/cache@v3
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

src/generalized.jl

@@ -60,6 +60,8 @@ pFq(α::NTuple{2, Any}, β::NTuple{1}, z; kwds...) = _₂F₁(α[1], α[2], β[1
 
 function pFq(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q}
     z = float(z)
+    n = terminating_order(α, β)
+    !isnothing(n) && return pFqmaclaurin(α, β, z; kwds...)
     if p ≤ q
         if real(z) > 0
             return pFqmaclaurin(α, β, z; kwds...)

src/specialfunctions.jl

@@ -153,6 +153,40 @@ end
 speciallogseries(x::Float64) = @clenshaw(5.0x, 1.0087391788544393911192, 1.220474262857857637288e-01, 8.7957928919918696061703e-03, 6.9050958578444820505037e-04, 5.7037120050065804396306e-05, 4.8731405131379353370205e-06, 4.2648797509486828820613e-07, 3.800372208946157617901e-08, 3.434168059359993493634e-09, 3.1381484326392473547608e-10, 2.8939845618385022798906e-11, 2.6892186934806386106143e-12, 2.5150879096374730760324e-13, 2.3652490233687788117887e-14, 2.2349973917002118259929e-15, 2.120769988408948118084e-16)
 speciallogseries(x::ComplexF64) = @evalpoly(x, 1.0000000000000000000000, 5.9999999999999999999966e-01, 4.2857142857142857142869e-01, 3.3333333333333333333347e-01, 2.7272727272727272727292e-01, 2.3076923076923076923072e-01, 1.9999999999999999999996e-01, 1.7647058823529411764702e-01, 1.5789473684210526315786e-01, 1.4285714285714285714283e-01, 1.3043478260869565217384e-01, 1.2000000000000000000000e-01, 1.1111111111111111111109e-01, 1.0344827586206896551722e-01, 9.6774193548387096774217e-02, 9.0909090909090909090938e-02, 8.5714285714285714285696e-02, 8.1081081081081081081064e-02, 7.6923076923076923076907e-02, 7.3170731707317073170688e-02)
 
+function terminating_order(α::Tuple, β::Tuple)
+    n = nothing
+    for αi in α
+        if αi ∈ ℤ && real(αi) ≤ 0
+            m = -round(Int, real(αi))
+            n = isnothing(n) ? m : min(n, m)
+        end
+    end
+    isnothing(n) && return nothing
+    for βi in β
+        if βi ∈ ℤ && real(βi) ≤ 0 && -round(Int, real(βi)) < n
+            return nothing
+        end
+    end
+    return n
+end
+
+function pFqterminating(α::NTuple{p, S}, β::NTuple{q, U}, z::V, n::Int) where {p, q, S, U, V}
+    T = promote_type(eltype(α), eltype(β), V)
+    term = one(T)
+    sum = one(T)
+    for k = 0:n-1
+        term *= z / (k + one(T))
+        for i = 1:p
+            term *= α[i] + k
+        end
+        for i = 1:q
+            term /= β[i] + k
+        end
+        sum += term
+    end
+    return sum
+end
+
 
 const libm = OpenLibm_jll.libopenlibm
 
@@ -537,6 +571,8 @@ function pFqmaclaurin(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where
 end
 
 function pFqmaclaurin(a::NTuple{p, S}, b::NTuple{q, U}, z::V; kmax::Int = KMAX) where {p, q, S, U, V}
+    n = terminating_order(a, b)
+    !isnothing(n) && return pFqterminating(a, b, z, n)
     T = promote_type(eltype(a), eltype(b), V)
     S₀, S₁, k = one(T), one(T)+prod(a)*z/prod(b), 1
     while k < kmax && errcheck(S₀, S₁, 8eps(real(T)))

test/confluent_tests.jl

@@ -0,0 +1,49 @@
+@testset "M" begin
+    @test M(-3, -3, 0.5) ≡ exp(0.5)
+    @test M(0, -1, 10) ≡ 1.0
+    @test M(1.2,  0.0, 0.1) ==  Inf # Mimick gamma( 0.0) =  Inf
+    @test M(1.2, -0.0, 0.1) == -Inf # and    gamma(-0.0) = -Inf
+    @test_throws DomainError M(1, -2, 0.5)
+    @test_throws DomainError M(-3, -2, 0.5)
+    @test M(-2, -3, 0.5) ≡ 1.375
+    @test M(0.5, 1.5, -1000) ≈ 0.028024956081989644 # From #46
+    @test M(1, 2, 0) == 1
+    @test M(1, 2, 0.25) == expm1(0.25)/0.25
+    for (S, T) in ((Float64, BigFloat),)
+        a = T(8.9)
+        b = T(0.5)
+        for x in T(-36):T(2):T(70)
+            @test M(S(a), S(b), S(x)) ≈ S(M(a, b, x)) # From #45
+        end
+        a = T(5)/6
+        b = T(1)/2
+        for x in T(-5):T(0.25):T(5)
+            @test M(S(a), S(b), -S(x)^2) ≈ S(M(a, b, -x^2)) # From #66
+        end
+        b = 1
+        x = T(1)/3
+        for a in S(1):S(0.5):S(7)
+            @test M(a, b, S(x)) ≈ S(M(a, b, x))
+        end
+    end
+end
+
+@testset "Confluent near-singular parameters" begin
+    for δ in (1e-8, 1e-10)
+        @test M(1.2, 1 + δ, 0.3) ≈ Float64(M(big(1.2), big(1 + δ), big(0.3)))
+        @test M(-2.5, -1 + δ, 0.2) ≈ Float64(M(big(-2.5), big(-1 + δ), big(0.2)))
+        @test _₂F₁(0.8, 0.7, 1.5 + δ, 0.9) ≈ Float64(_₂F₁(big(0.8), big(0.7), big(1.5 + δ), big(0.9)))
+    end
+end
+
+@testset "U" begin
+    @test U(1, 1, 1.f0) ≈ 0.5963473623231942 # the Euler series
+    @test U(1, 1, 1) == 0.5963473623231942
+    for (S, T) in ((Float64, BigFloat),)
+        b = 0
+        x = T(1)/3
+        for a in S(1):S(0.5):S(7)
+            @test U(a, b, S(x)) ≈ S(U(a, b, x)) # From #55
+        end
+    end
+end

test/regression_tests.jl

@@ -0,0 +1,85 @@
+@testset "Exact identities" begin
+    z = 0.3
+    @test M(2.5, 2.5, z) ≈ exp(z)
+    @test pFq((3,), (), z) ≈ (1 - z)^(-3)
+    @test pFq([1.0, 2.0], [4.0], 0.25) == pFq((1.0, 2.0), (4.0,), 0.25)
+    @test _₂F₁(0, 2, 3, 0.7) == 1.0
+    @test _₂F₁(2, 3, 2, 0.25) == pFq((3,), (), 0.25)
+    @test _₂F₁(2, 3, 3, 0.25) == pFq((2,), (), 0.25)
+end
+
+@testset "Singular parameters" begin
+    @test_throws DomainError M(1.2, -1, 0.1)
+    @test isinf(_₂F₁(1, 2, 2, 1.0))
+    @test isnan(_₂F₁(1, 2 + im, 3, 1))
+end
+
+@testset "Polyalgorithm boundaries" begin
+    α, β = (1.0, 2.0), (3.0,)
+    @test pFq(α, β, 0.72) ≈ pFq(Tuple(α), Tuple(β), 0.72)
+    @test pFq(α, β, 0.73) ≈ pFq(Tuple(α), Tuple(β), 0.73)
+    @test pFq((1.0,), (2.0,), 0.0) == 1.0
+    @test pFq((1.0, 2.0), (3.0,), 0.0) == 1.0
+end
+
+@testset "₂F₁ special cases" begin
+    @test _₂F₁(0.3, -0.3, 0.5, 0.2) ≈ HypergeometricFunctions.cosnasinsqrt(-0.6, 0.2)
+    @test _₂F₁(0.5, 0.5, 1.5, 0.2) ≈ HypergeometricFunctions.sqrtasinsqrt(0.2)
+    @test _₂F₁(0.5, 1.0, 1.5, 0.2) ≈ HypergeometricFunctions.sqrtatanhsqrt(0.2)
+    @test _₂F₁(1.0, 1.0, 2.0, 0.2) ≈ HypergeometricFunctions.log1pover(-0.2)
+    @test _₂F₁(2.0, 2.0, 4.0, 0.2) ≈ HypergeometricFunctions.logandpoly(0.2)
+    @test _₂F₁(1.0, 1.5, 2.5, 0.2) ≈ HypergeometricFunctions.speciallog(0.2)
+end
+
+@testset "Branch cut consistency" begin
+    for ϵ in (1e-8, 1e-10)
+        z₊ = -0.3 + im * ϵ
+        z₋ = -0.3 - im * ϵ
+        @test pFq((0.5,), (), z₊) == conj(pFq((0.5,), (), z₋))
+
+        z₊ = 1.2 + im * ϵ
+        z₋ = 1.2 - im * ϵ
+        @test M(0.75, 1.25, z₊) == conj(M(0.75, 1.25, z₋))
+
+        z₊ = 0.9 + im * ϵ
+        z₋ = 0.9 - im * ϵ
+        @test _₂F₁(0.8, 0.7, 1.5, z₊) == conj(_₂F₁(0.8, 0.7, 1.5, z₋))
+    end
+end
+
+@testset "Issue regressions" begin
+    function terminating_reference(α, β, z)
+        n = minimum(-round(Int, real(ai)) for ai in α if isinteger(ai) && real(ai) ≤ 0)
+        term = one(promote_type(eltype(α), eltype(β), typeof(z)))
+        sum = term
+        for k = 0:n-1
+            term *= z / (k + 1)
+            for ai in α
+                term *= ai + k
+            end
+            for bi in β
+                term /= bi + k
+            end
+            sum += term
+        end
+        return sum
+    end
+
+    @test pFq((-4, -3, 151), (2, -153), -1.0) ≈ terminating_reference((-4, -3, 151), (2, -153), -1.0)
+
+    for (a, b) in (
+        ((-9.0, 16.5, 12.5), (29.0, -32.5)),
+        ((-13.0, 4.5, 3.5), (8.0, -17.5)),
+        ((-21.0, 3.5, 2.5), (6.0, -23.5)),
+        ((-29.0, 3.5, 2.5), (6.0, -32.5)),
+    )
+        @test pFq(a, b, 1.0) ≈ terminating_reference(a, b, 1.0) rtol = 1e-12
+    end
+
+    @test _₂F₁(-1, 0, -1, 2.0) == 1.0
+    @test _₂F₁(-1, 2, -1, 2.0) == pFq((2,), (), 2.0)
+end
+
+@testset "Warning symbols" begin
+    @test pFq2string(Val(86420), Val(97531)) == "₈₆₄₂₀F₉₇₅₃₁"
+end