close #66

Author: MikaelSlevinsky

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"

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

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...)

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.

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...)

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