diff --git a/Project.toml b/Project.toml index dc2af52..c58c18f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "HypergeometricFunctions" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.17" +version = "0.3.18" [deps] DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" diff --git a/README.md b/README.md index 66061f7..011820f 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ julia> pFq((1, ), (2, ), 0.01) # ≡ expm1(0.01)/0.01 1.0050167084168058 julia> pFq((1/3, ), (2/3, ), -1000) # ₁F₁ -0.050558053946448994 +0.050558053946448855 julia> pFq((1, 2), (4, ), 1) # a well-poised ₂F₁ 2.9999999999999996 diff --git a/src/confluent.jl b/src/confluent.jl index 40a1497..dd2f9c4 100644 --- a/src/confluent.jl +++ b/src/confluent.jl @@ -1,7 +1,7 @@ # The references to special cases are to NIST's DLMF. """ -Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a; b; z)`. +Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a, b; z)`. """ function _₁F₁(a, b, z; kwds...) z = float(z) @@ -24,30 +24,20 @@ function _₁F₁(a, b, z; kwds...) end function _₁F₁general(a, b, z; kwds...) - if abs(z) > 10 # TODO: check this algorithmic parameter cutoff, determine when to use rational algorithms - if isreal(z) - if real(z) ≥ 0 # 13.7.1 - return gamma(b)/gamma(a)*exp(z)*z^(a-b)*pFq((1-a, b-a), (), inv(z); kwds...) - else # 13.7.1 + 13.2.39 - return gamma(b)/gamma(b-a)*(-z)^-a*pFq((a, a-b+1), (), -inv(z); kwds...) - end - else # 13.7.2 - return gamma(b)/gamma(a)*exp(z)*z^(a-b)*pFq((1-a, b-a), (), inv(z); kwds...) + gamma(b)/gamma(b-a)*(-z)^-a*pFq((a, a-b+1), (), -inv(z); kwds...) - end - elseif real(z) ≥ 0 + if real(z) > 0 return _₁F₁maclaurin(a, b, z; kwds...) - else # 13.2.39 - return exp(z)*_₁F₁(b-a, b, -z; kwds...) + else + return pFqweniger((a, ), (b, ), z; kwds...) end end """ -Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a; b; z)`. +Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a, b; z)`. """ const M = _₁F₁ """ -Compute Tricomi's confluent hypergeometric function `U(a, b, z) ∼ z⁻ᵃ ₂F₀([a, a-b+1]; []; -z⁻¹)`. +Compute Tricomi's confluent hypergeometric function `U(a, b, z) ∼ z⁻ᵃ ₂F₀((a, a-b+1), (); -z⁻¹)`. """ function U(a, b, z; kwds...) return z^-a*pFq((a, a-b+1), (), -inv(z); kwds...) diff --git a/src/gauss.jl b/src/gauss.jl index 6c96e8b..ef55888 100644 --- a/src/gauss.jl +++ b/src/gauss.jl @@ -1,7 +1,7 @@ # The references to special cases are to Table of Integrals, Series, and Products, § 9.121, followed by NIST's DLMF. """ -Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)`. +Compute the Gauss hypergeometric function `₂F₁(a, b, c; z)`. """ function _₂F₁(a, b, c, z; method::Symbol = :general, kwds...) z = float(z) @@ -56,7 +56,7 @@ function _₂F₁(a, b, c, z; method::Symbol = :general, kwds...) end """ -Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with positive parameters a, b, and c and argument 0 ≤ z ≤ 1. Useful for statisticians. +Compute the Gauss hypergeometric function `₂F₁(a, b, c; z)` with positive parameters a, b, and c and argument 0 ≤ z ≤ 1. Useful for statisticians. """ function _₂F₁positive(a, b, c, z; kwds...) @assert a > 0 && b > 0 && c > 0 && 0 ≤ z ≤ 1 @@ -68,7 +68,7 @@ function _₂F₁positive(a, b, c, z; kwds...) end """ -Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with general parameters a, b, and c. +Compute the Gauss hypergeometric function `₂F₁(a, b, c; z)` with general parameters a, b, and c. This polyalgorithm is designed based on the paper N. Michel and M. V. Stoitsov, Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl–Teller–Ginocchio potential wave functions, Comp. Phys. Commun., 178:535–551, 2008. @@ -98,7 +98,7 @@ function _₂F₁general(a, b, c, z; kwds...) end """ -Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with general parameters a, b, and c. +Compute the Gauss hypergeometric function `₂F₁(a, b, c; z)` with general parameters a, b, and c. This polyalgorithm is designed based on the review J. W. Pearson, S. Olver and M. A. Porter, Numerical methos for the computation of the confluent and Gauss hypergeometric functions, arXiv:1407.7786, 2014. diff --git a/src/generalized.jl b/src/generalized.jl index cf550ae..b77c3d9 100644 --- a/src/generalized.jl +++ b/src/generalized.jl @@ -15,7 +15,7 @@ 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) if p ≤ q - if real(z) ≥ 0 + if real(z) > 0 return pFqmaclaurin(α, β, z; kwds...) else return pFqweniger(α, β, z; kwds...) @@ -45,8 +45,11 @@ function pFqcontinuedfraction(α::AbstractVector{S}, β::AbstractVector{U}, z::V return 1 + z * prod(α) / prod(β) / (1 + K) end +""" +Compute the generalized hypergeometric function `₃F₂(a₁, 1, 1, b₁, 2; z)`. +""" +_₃F₂(a₁, b₁, z; kwds...) = _₃F₂(a₁, 1, 1, b₁, 2, z; kwds...) """ Compute the generalized hypergeometric function `₃F₂(a₁, a₂, a₃, b₁, b₂; z)`. """ _₃F₂(a₁, a₂, a₃, b₁, b₂, z; kwds...) = pFq((a₁, a₂, a₃), (b₁, b₂), z; kwds...) -_₃F₂(a₁, b₁, z; kwds...) = _₃F₂(a₁, 1, 1, b₁, 2, z; kwds...) diff --git a/test/runtests.jl b/test/runtests.jl index 9687069..5b9c782 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -480,23 +480,32 @@ end @test M(-2, -3, 0.5) ≡ 1.375 @test M(0.5, 1.5, -1000) ≈ 0.028024956081989644 # From #46 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 - z = T(1)/3 - x = S(z) + x = T(1)/3 for a in S(1):S(0.5):S(7) - @test M(a, b, x) ≈ S(M(a, b, z)) + @test M(a, b, S(x)) ≈ S(M(a, b, x)) end end end @testset "U" begin - # From #55 + @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 - z = T(1)/3 - x = S(z) + x = T(1)/3 for a in S(1):S(0.5):S(7) - @test U(a, b, x) ≈ S(U(a, b, z)) + @test U(a, b, S(x)) ≈ S(U(a, b, x)) # From #55 end end end