let method=:general2 be valid

MikaelSlevinsky · MikaelSlevinsky · commit 5813051d4c5f · 2023-07-22T22:59:37.000-06:00
diff --git a/src/gauss.jl b/src/gauss.jl
@@ -50,6 +50,8 @@ function _₂F₁(a, b, c, z; method::Symbol = :general, kwds...)
     end
     if method == :positive
         return _₂F₁positive(a, b, c, z; kwds...)
+    elseif method == :general2
+        return _₂F₁general2(a, b, c, z; kwds...)
     else#if method == :general
         return _₂F₁general(a, b, c, z; kwds...) # catch-all
     end
@@ -103,55 +105,55 @@ This polyalgorithm is designed based on the review
 
 > J. W. Pearson, S. Olver and M. A. Porter, [Numerical methods for the computation of the confluent and Gauss hypergeometric functions](https://doi.org/10.1007/s11075-016-0173-0), *Numer. Algor.*, **74**:821–866, 2017.
 """
-function _₂F₁general2(a, b, c, z)
+function _₂F₁general2(a, b, c, z; kwds...)
     T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
     if z == 1
-        return _₂F₁argument_unity(a, b, c, z)
+        return _₂F₁argument_unity(a, b, c, z; kwds...)
     elseif real(b) < real(a)
-        return _₂F₁general2(b, a, c, z)
+        return _₂F₁general2(b, a, c, z; kwds...)
     elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
-        return _₂F₁maclaurin(a, b, c, z)
+        return _₂F₁maclaurin(a, b, c, z; kwds...)
     elseif !isalmostwellpoised(a, b, c)
-        return exp((c-a-b)*log1p(-z))*_₂F₁general2(c-a, c-b, c, z)
+        return exp((c-a-b)*log1p(-z))*_₂F₁general2(c-a, c-b, c, z; kwds...)
     elseif abs(z / (z - 1)) ≤ ρ && absarg(1 - z) < convert(real(T), π) # 15.8.1
         w = z/(z-1)
-        return _₂F₁maclaurin(a, c-b, c, w)*exp(-a*log1p(-z))
+        return _₂F₁maclaurin(a, c-b, c, w; kwds...)*exp(-a*log1p(-z))
     elseif abs(inv(z)) ≤ ρ && absarg(-z) < convert(real(T), π)
         w = inv(z)
         if isapprox(a, b) # 15.8.8
-            return gamma(c)/gamma(a)/gamma(c-a)*(-w)^a*_₂F₁logsumalt(a, c-a, z, w)
+            return gamma(c)/gamma(a)/gamma(c-a)*(-w)^a*_₂F₁logsumalt(a, c-a, z, w; kwds...)
         elseif a-b ∉ ℤ # 15.8.2
-            return gamma(c)*((-w)^a*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, a-c+1, a-b+1, w)+(-w)^b*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, b-c+1, b-a+1, w))
+            return gamma(c)*((-w)^a*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, a-c+1, a-b+1, w; kwds...)+(-w)^b*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, b-c+1, b-a+1, w; kwds...))
         end # TODO: full 15.8.8
     elseif abs(inv(1-z)) ≤ ρ && absarg(-z) < convert(real(T), π)
         w = inv(1-z)
         if isapprox(a, b) # 15.8.9
-            return gamma(c)*exp(-a*log1p(-z))/gamma(a)/gamma(c-b)*_₂F₁logsum(a, c-a, z, w, 1)
+            return gamma(c)*exp(-a*log1p(-z))/gamma(a)/gamma(c-b)*_₂F₁logsum(a, c-a, z, w, 1; kwds...)
         elseif a-b ∉ ℤ # 15.8.3
-            return gamma(c)*(exp(-a*log1p(-z))*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, c-b, a-b+1, w)+exp(-b*log1p(-z))*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, c-a, b-a+1, w))
+            return gamma(c)*(exp(-a*log1p(-z))*gamma(b-a)/gamma(b)/gamma(c-a)*_₂F₁maclaurin(a, c-b, a-b+1, w; kwds...)+exp(-b*log1p(-z))*gamma(a-b)/gamma(a)/gamma(c-b)*_₂F₁maclaurin(b, c-a, b-a+1, w; kwds...))
         end # TODO: full 15.8.9
     elseif abs(1-z) ≤ ρ && absarg(z) < convert(real(T), π) && absarg(1-z) < convert(real(T), π)
         w = 1-z
         if isapprox(c, a + b) # 15.8.10
-            return gamma(c)/gamma(a)/gamma(b)*_₂F₁logsum(a, b, z, w, -1)
+            return gamma(c)/gamma(a)/gamma(b)*_₂F₁logsum(a, b, z, w, -1; kwds...)
         elseif c - a - b ∉ ℤ # 15.8.4
-            return gamma(c)*(gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, b, a+b-c+1, w)+exp((c-a-b)*log1p(-z))*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, c-b, c-a-b+1, w))
+            return gamma(c)*(gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, b, a+b-c+1, w; kwds...)+exp((c-a-b)*log1p(-z))*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, c-b, c-a-b+1, w; kwds...))
         end # TODO: full 15.8.10
     elseif abs(1-inv(z)) ≤ ρ && absarg(z) < convert(real(T), π) && absarg(1-z) < convert(real(T), π)
         w = 1-inv(z)
         if isapprox(c, a + b) # 15.8.11
-            return gamma(c)*z^(-a)/gamma(a)*_₂F₁logsumalt(a, b, z, w)
+            return gamma(c)*z^(-a)/gamma(a)*_₂F₁logsumalt(a, b, z, w; kwds...)
         elseif c - a - b ∉ ℤ # 15.8.5
-            return gamma(c)*(z^(-a)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, a-c+1, a+b-c+1, w)+z^(a-c)*(1-z)^(c-a-b)*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, 1-a, c-a-b+1, w))
+            return gamma(c)*(z^(-a)*gamma(c-a-b)/gamma(c-a)/gamma(c-b)*_₂F₁maclaurin(a, a-c+1, a+b-c+1, w; kwds...)+z^(a-c)*(1-z)^(c-a-b)*gamma(a+b-c)/gamma(a)/gamma(b)*_₂F₁maclaurin(c-a, 1-a, c-a-b+1, w; kwds...))
         end # TODO: full 15.8.11
     elseif abs(z-0.5) > 0.5
         if isapprox(a, b) && !isapprox(c, a+0.5)
-            return gamma(c)/gamma(a)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuationalt(a, c, 0.5, z)
+            return gamma(c)/gamma(a)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuationalt(a, c, 0.5, z; kwds...)
         elseif a-b ∉ ℤ
-            return gamma(c)*(gamma(b-a)/gamma(b)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuation(a, a+b, c, 0.5, z) + gamma(a-b)/gamma(a)/gamma(c-b)*(0.5-z)^(-b)*_₂F₁continuation(b, a+b, c, 0.5, z))
+            return gamma(c)*(gamma(b-a)/gamma(b)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuation(a, a+b, c, 0.5, z; kwds...) + gamma(a-b)/gamma(a)/gamma(c-b)*(0.5-z)^(-b)*_₂F₁continuation(b, a+b, c, 0.5, z; kwds...))
         end
     end
-    return pFqweniger((a, b), (c, ), z)
+    return pFqweniger((a, b), (c, ), z; kwds...) # _₂F₁taylor(a, b, c, z; kwds...)
 end
 
 # Special case of (-x)^a*_₂F₁ to handle LogNumber correctly in RiemannHilbert.jl
diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl
@@ -471,7 +471,7 @@ function _₂F₁continuationalt(a::Number, c::Number, z₀::Number, z::Number)
     izz₀ = inv(z-z₀)
     e0, e1 = one(T), (a+one(T))*(one(T)-2z₀)+(2a+one(T))*z₀-c
     f0, f1 = zero(T), one(T)-2z₀
-    cⱼ = log(z₀-z)+2digamma(one(T))-digamma(a)-digamma(c-a)
+    cⱼ = log(z₀-z)+2digamma(one(real(T)))-digamma(a)-digamma(c-a)
     S₀ = cⱼ
     cⱼ += 2/one(T)-one(T)/a
     C = a*izz₀
@@ -489,7 +489,7 @@ end
 
 function _₂F₁logsum(a::Number, b::Number, z::Number, w::Number, s::Int)
     T = promote_type(typeof(a), typeof(b), typeof(z), typeof(w))
-    cⱼ = 2digamma(one(T))-digamma(a)-digamma(b)+s*log1p(-z)
+    cⱼ = 2digamma(one(real(T)))-digamma(a)-digamma(b)+s*log1p(-z)
     C, S, j = one(T), cⱼ, 0
     while abs(C) > 8abs(S)*eps(real(T)) || j ≤ 1
         C *= (a+j)/(j+1)^2*(b+j)*w
@@ -502,7 +502,7 @@ end
 
 function _₂F₁logsumalt(a::Number, b::Number, z::Number, w::Number)
     T = promote_type(typeof(a), typeof(b), typeof(z), typeof(w))
-    d, cⱼ = one(T)-b, 2digamma(one(T))-digamma(a)-digamma(b)-log(-w)
+    d, cⱼ = one(T)-b, 2digamma(one(real(T)))-digamma(a)-digamma(b)-log(-w)
     C, S, j = one(T), cⱼ, 0
     while abs(C) > 8abs(S)*eps(real(T)) || j ≤ 1
         C *= (a+j)/(j+1)^2*(d+j)*w
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -1,8 +1,8 @@
 using HypergeometricFunctions, SpecialFunctions, Test
 import LinearAlgebra: norm
 import HypergeometricFunctions: iswellpoised, isalmostwellpoised, M, U,
-                                pochhammer,
-                                _₂F₁general, pFqdrummond, pFqweniger, pFq2string
+                                pochhammer, _₂F₁general, _₂F₁general2,
+                                pFqdrummond, pFqweniger, pFq2string
 
 const rtol = 1.0e-3
 const NumberType = Float64
@@ -24,9 +24,9 @@ const NumberType = Float64
 end
 
 @testset "Hypergeometric Functions" begin
-    @testset "_₂F₁ vs _₂F₁general" begin
+    @testset "_₂F₁ vs _₂F₁general vs _₂F₁general2" begin
         e = exp(1.0)
-        regression_max_accumulated_error = 8.0e-13 # discovered by running the test
+        regression_max_accumulated_error = 2^12*eps() # discovered by running the test
         for z in (.9rand(Float64, 10), 10rand(ComplexF64, 10))
             j = 1
             for (a,b,c) in ((√2/e, 1.3, 1.3), (1.2, √3, 1.2), (-0.4, 0.4, 0.5),
@@ -46,10 +46,27 @@ end
                 @test error_accum < regression_max_accumulated_error
                 j += 1
             end
+            for (a,b,c) in ((√2/e, 1.3, 1.3), (1.2, √3, 1.2), (-0.4, 0.4, 0.5),
+                            (-0.3, 1.3, 0.5), (0.35, 0.85, 0.5), (0.5, 1.0, 1.5),
+                            (0.3, 0.7, 1.5), (0.7, 1.3, 1.5), (0.35, 0.85, 1.5),
+                            (3.0, 1.0, 2.0), (-2.0, 1.0, 2.0), (-3.0, 1.0, 2.0),
+                            (1.0, -4.0, 2.0), (1.0, 1.5, 2.5))
+                error_accum = 0.0
+                for zi in z
+                    twoFone = _₂F₁(a,b,c,zi)
+                    aa,bb,cc = big(a),big(b),big(c)
+                    twoFonegeneral2 = convert(Complex{Float64},_₂F₁general2(aa, bb, cc, big(zi)))
+                    norm(twoFone / twoFonegeneral2 - 1) > sqrt(eps()) && println("This is ₂F₁($a,$b;$c;zi) - ₂F₁general2($a,$b;$c;zi): ", norm(twoFone / twoFonegeneral2 - 1), "   ", twoFone, "   ", twoFonegeneral2, "   ", isfinite(twoFone), "   ", isfinite(twoFonegeneral2), " this is zi: ", zi)
+                    error_accum += Float64(norm(twoFone / twoFonegeneral2 - 1))
+                end
+                @test error_accum < regression_max_accumulated_error
+                j += 1
+            end
         end
-        @testset "Test that _₂F₁ is inferred for Float32 arguments" begin
-            @test @inferred(_₂F₁(0.3f0, 0.7f0, 1.3f0, 0.1f0)) ≈ Float32(_₂F₁(0.3, 0.7, 1.3, 0.1))
-        end
+    end
+
+    @testset "Test that _₂F₁ is inferred for Float32 arguments" begin
+        @test @inferred(_₂F₁(0.3f0, 0.7f0, 1.3f0, 0.1f0)) ≈ Float32(_₂F₁(0.3, 0.7, 1.3, 0.1))
     end
 
     @testset "method = positive" begin
