Consistent evaluation of Tricomi U(a,b,z)

src/confluent.jl

@@ -40,5 +40,61 @@ const M = _₁F₁
 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...)
+    z = float(z)
+    if z isa Real && z < 0
+        z = complex(z)
+    end
+
+    if isequal(a, 1) && isequal(b, 1)
+        return _U11(z)
+    elseif -a ∈ ℕ₀
+        return _U_polynomial(-round(Int, real(a)), b, z)
+    elseif b ∈ ℤ
+        n = round(Int, real(b))
+        if n > 0
+            return _U_integer_b(a, n - 1, z)
+        else
+            return exp((1 - b) * log(z)) * _U_integer_b(a - b + 1, 1 - n, z)
+        end
+    else
+        return gamma(1 - b) * ogamma(a - b + 1) * _₁F₁(a, b, z; kwds...) +
+               gamma(b - 1) * ogamma(a) * exp((1 - b) * log(z)) * _₁F₁(a - b + 1, 2 - b, z; kwds...)
+    end
+end
+
+_U11(z::Float32) = Float32(_U11(Float64(z)))
+_U11(z::ComplexF32) = ComplexF32(_U11(ComplexF64(z)))
+_U11(z) = exp(z) * expint(z)
+
+function _U_polynomial(m::Int, b, z)
+    T = promote_type(typeof(b), typeof(z))
+    ret = zero(T)
+    for s = 0:m
+        ret += (-1)^s * binomial(m, s) * pochhammer(b + s, m - s) * z^s
+    end
+    return ret
+end
+
+function _U_integer_b(a, n::Int, z)
+    T = promote_type(typeof(a), typeof(z))
+    logz = log(z)
+    sum1 = zero(T)
+    term = one(T)
+    tol = 10eps(real(T))
+    for k = 0:KMAX
+        coeff = logz + digamma(a + k) - digamma(k + 1) - digamma(n + k + 1)
+        addend = term * coeff
+        sum1 += addend
+        if k > 0 && !errcheck(sum1, sum1 - addend, tol)
+            break
+        end
+        term *= (a + k) * z / ((n + k + 1) * (k + 1))
+    end
+
+    sum2 = zero(T)
+    for k = 1:n
+        sum2 += factorial(k - 1) * pochhammer(1 - a + k, n - k) * z^(-k) / factorial(n - k)
+    end
+
+    return (-1)^(n + 1) * ogamma(a - n) / factorial(n) * sum1 + ogamma(a) * sum2
 end

test/runtests.jl

@@ -601,14 +601,26 @@ end
 
 @testset "U" begin
     @test U(1, 1, 1.f0) ≈ 0.5963473623231942 # the Euler series
-    @test U(1, 1, 1) == 0.5963473623231942
+    @test U(1, 1, 1) ≈ 0.5963473623231942
+    @test U(1, 1, 1 + 0.2im) ≈ exp(1 + 0.2im) * expint(1 + 0.2im)
+    @test U(1, 1, ComplexF32(1 + 0.2im)) ≈ ComplexF32(exp(ComplexF64(1 + 0.2im)) * expint(ComplexF64(1 + 0.2im)))
+    @test U(-2, 3/2, 0.3) ≈ 0.3^2 - 2 * (3/2 + 1) * 0.3 + (3/2) * (3/2 + 1)
+    @test U(1, 2, 0.3) ≈ inv(0.3)
+    @test U(3/2, 0, 0.3) ≈ 0.3 * U(5/2, 2, 0.3)
     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
+
+    # Regression tests from #63, which is related to the negative-real-axis
+    # branch behavior discussed in #61.
+    @test U(0.5, 0.5, Complex(-2.5)) ≈ 0.145492 - 0.810467im atol=1e-6 rtol=1e-6
+    @test U(-0.5, -0.5, Complex(-2.5)) ≈ 0.0727459 + 1.17591im atol=1e-5 rtol=1e-5
+    @test U(-0.5, -0.5, -2.5 - 1.5im) ≈ 0.595921 - 1.34976im atol=1e-5 rtol=1e-5
+    @test U(0.5, 0.5, -2.5) ≈ U(0.5, 0.5, Complex(-2.5))
 end
 
 @testset "Warning symbols" begin