add Tricomi's U through rational approximation of its asymptotic expa…

…nsion

Project.toml

@@ -1,6 +1,6 @@
 name = "HypergeometricFunctions"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
-version = "0.3.7"
+version = "0.3.8"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"

README.md

@@ -4,4 +4,4 @@
 
 A Julia package for calculating hypergeometric functions
 
-This package implements the generalized hypergeometric function `pFq([a1,…,am], [b1,…,bn], z)`. In particular, the Gauss hypergeometric function is available as `_₂F₁(a, b, c, z)`, Kummer's confluent hypergeometric function is available as `_₁F₁(a, b, z)` and `HypergeometricFunctions.M(a, b, z)`, as well as `_₃F₂([a1, a2, a3], [b1, b2], z)`.
+This package implements the generalized hypergeometric function `pFq([a1,…,am], [b1,…,bn], z)`. In particular, the Gauss hypergeometric function is available as `_₂F₁(a, b, c, z)`, confluent hypergeometric function is available as `_₁F₁(a, b, z) ≡ HypergeometricFunctions.M(a, b, z)` and `HypergeometricFunctions.U(a, b, z)`, as well as `_₃F₂([a1, a2, a3], [b1, b2], z)`.

src/HypergeometricFunctions.jl

@@ -11,7 +11,7 @@ export _₁F₁, _₂F₁, _₃F₂, pFq
 
 include("specialfunctions.jl")
 include("gauss.jl")
-include("kummer.jl")
+include("confluent.jl")
 include("generalized.jl")
 include("drummond.jl")
 include("weniger.jl")

src/confluent.jl

@@ -0,0 +1,22 @@
+"""
+Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a; b; z)`.
+"""
+function _₁F₁(a, b, z)
+    if real(z) ≥ 0
+        return _₁F₁maclaurin(a, b, z)
+    else
+        return exp(z)*_₁F₁(b-a, b, -z)
+    end
+end
+
+"""
+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⁻¹)`.
+"""
+function U(a, b, z)
+    return z^-a*pFq([a, a-b+1], typeof(z)[], -inv(z))
+end

src/generalized.jl

@@ -14,8 +14,18 @@ function pFq(α::AbstractVector, β::AbstractVector, z; kwds...)
         return _₁F₁(α[1], β[1], float(z))
     elseif length(α) == 2 && length(β) == 1
         return _₂F₁(α[1], α[2], β[1], float(z))
-    elseif abs(z) ≤ ρ
-        return pFqmaclaurin(α, β, float(z); kwds...)
+    elseif length(α) ≤ length(β)
+        if real(z) ≥ 0
+            return pFqmaclaurin(α, β, float(z); kwds...)
+        else
+            return pFqweniger(α, β, float(z); kwds...)
+        end
+    elseif length(α) == length(β) + 1
+        if abs(z) ≤ ρ
+            return pFqmaclaurin(α, β, float(z); kwds...)
+        else
+            return pFqweniger(α, β, float(z); kwds...)
+        end
     else
         return pFqweniger(α, β, float(z); kwds...)
     end