diff --git a/.codecov.yml b/.codecov.yml
deleted file mode 100644
index 8b13789..0000000
--- a/.codecov.yml
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/Project.toml b/Project.toml
index b6ead10..ac8e51e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,21 +1,14 @@
 name = "HypergeometricFunctions"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
-version = "0.2.3"
+version = "0.3"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-DualNumbers = "0.5, 0.6"
-IntervalArithmetic = "0.15, 0.16"
+DualNumbers = "0.5, 0.6, 0.7"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10"
 julia = "1"
-
-[extras]
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test", "IntervalArithmetic"]
diff --git a/README.md b/README.md
index 06f9442..69c2813 100644
--- a/README.md
+++ b/README.md
@@ -5,4 +5,4 @@
 # HypergeometricFunctions.jl
 A Julia package for calculating hypergeometric functions
 
-This package implements some hypergeometric functions. In particular, the Gauss hypergeometric function is available as `_₂F₁(a,b,c,z)`, and also `_₃F₂([a1, a2, a3], [b1,b2], z)` and more general `mFn([a1,…,am], [b1,…,bn], 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)`, and also `_₃F₂([a1, a2, a3], [b1, b2], z)`.
diff --git a/src/HypergeometricFunctions.jl b/src/HypergeometricFunctions.jl
index d71c8a2..a180a1f 100644
--- a/src/HypergeometricFunctions.jl
+++ b/src/HypergeometricFunctions.jl
@@ -1,16 +1,17 @@
 """
 This module calculates (generalized) hypergeometric functions:
 
-    mFn(a;b;z) = Σ_{k=0}^∞ (a_1,k) ⋯ (a_m,k) / (b_1,k) ⋯ (b_n,k) zᵏ/k!
+    pFq(α, β; z) = Σ_{k=0}^∞ (α_1)ₖ ⋯ (α_p)ₖ / (β_1)ₖ ⋯ (β_q)ₖ zᵏ/k!
 """
 module HypergeometricFunctions
 
 using DualNumbers, LinearAlgebra, SpecialFunctions
 
-export _₂F₁, _₃F₂, mFn
+export _₂F₁, _₃F₂, pFq
 
-include("gauss.jl")
 include("specialfunctions.jl")
+include("gauss.jl")
+include("generalized.jl")
 include("drummond.jl")
 include("weniger.jl")
 
diff --git a/src/drummond.jl b/src/drummond.jl
index a5bd3ae..1dd09ea 100644
--- a/src/drummond.jl
+++ b/src/drummond.jl
@@ -2,7 +2,7 @@
 # using Drummond's sequence transformation
 
 # ₀F₀(;z)
-function drummond0F0(z::T) where T
+function drummond0F0(z::T; kmax::Int = 10_000) where T
     if norm(z) < eps(real(T))
         return one(T)
     end
@@ -18,7 +18,7 @@ function drummond0F0(z::T) where T
     Dhi, Dlo = ((k+2)*ζ-1)*Dhi + k*ζ*Dlo, Dhi
     Thi, Tlo = Nhi/Dhi, Thi
     k += 1
-    while abs(Thi-Tlo) > 10*abs(Thi)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tlo, Thi, 10eps(real(T)))
         Nhi, Nlo = ((k+2)*ζ-1)*Nhi + k*ζ*Nlo, Nhi
         Dhi, Dlo = ((k+2)*ζ-1)*Dhi + k*ζ*Dlo, Dhi
         Thi, Tlo = Nhi/Dhi, Thi
@@ -28,9 +28,9 @@ function drummond0F0(z::T) where T
 end
 
 # ₁F₀(α;z)
-function drummond1F0(α::T1, z::T2) where {T1, T2}
+function drummond1F0(α::T1, z::T2; kmax::Int = 10_000) where {T1, T2}
     T = promote_type(T1, T2)
-    absα = T(abs(α))
+    absα = abs(T(α))
     if norm(z) < eps(real(T)) || norm(α) < eps(absα)
         return one(T)
     end
@@ -56,7 +56,7 @@ function drummond1F0(α::T1, z::T2) where {T1, T2}
     Nhi /= α+k+1
     Dhi /= α+k+1
     k += 1
-    while abs(Thi-Tlo) > 10*abs(Thi)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tlo, Thi, 10eps(real(T)))
         Nhi, Nlo = ((k+2)*ζ-(α+2k+1))*Nhi + k*(ζ-1)*Nlo, Nhi
         Dhi, Dlo = ((k+2)*ζ-(α+2k+1))*Dhi + k*(ζ-1)*Dlo, Dhi
         Thi, Tlo = Nhi/Dhi, Thi
@@ -71,7 +71,7 @@ function drummond1F0(α::T1, z::T2) where {T1, T2}
 end
 
 # ₀F₁(β;z)
-function drummond0F1(β::T1, z::T2) where {T1, T2}
+function drummond0F1(β::T1, z::T2; kmax::Int = 10_000) where {T1, T2}
     T = promote_type(T1, T2)
     if norm(z) < eps(real(T))
         return one(T)
@@ -91,7 +91,7 @@ function drummond0F1(β::T1, z::T2) where {T1, T2}
     Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-1)*Dhi + k*(β+2k+2)*ζ*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
     Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
     k += 1
-    while abs(Thi-Tmid) > 10*abs(Thi)*eps(real(T)) && abs(Tmid-Tlo) > 10*abs(Tmid)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
         Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-1)*Nhi + k*(β+2k+2)*ζ*Nmid + k*(k-1)*ζ*Nlo, Nhi, Nmid
         Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-1)*Dhi + k*(β+2k+2)*ζ*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
         Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
@@ -101,10 +101,10 @@ function drummond0F1(β::T1, z::T2) where {T1, T2}
 end
 
 # ₂F₀(α,β;z)
-function drummond2F0(α::T1, β::T2, z::T3) where {T1, T2, T3}
+function drummond2F0(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
     T = promote_type(T1, T2, T3)
-    absα = T(abs(α))
-    absβ = T(abs(β))
+    absα = abs(T(α))
+    absβ = abs(T(β))
     if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ)
         return one(T)
     end
@@ -129,7 +129,7 @@ function drummond2F0(α::T1, β::T2, z::T3) where {T1, T2, T3}
     Nhi /= (α+2)*(β+2)
     Dhi /= (α+2)*(β+2)
     k = 2
-    while abs(Thi-Tmid) > 10*abs(Thi)*eps(real(T)) && abs(Tmid-Tlo) > 10*abs(Tmid)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
         Nhi, Nmid, Nlo = ((k+2)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Nhi - k*(α+β+3k-ζ)*Nmid - k*(k-1)*Nlo, Nhi, Nmid
         Dhi, Dmid, Dlo = ((k+2)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Dhi - k*(α+β+3k-ζ)*Dmid - k*(k-1)*Dlo, Dhi, Dmid
         Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
@@ -144,9 +144,9 @@ function drummond2F0(α::T1, β::T2, z::T3) where {T1, T2, T3}
 end
 
 # ₁F₁(α,β;z)
-function drummond1F1(α::T1, β::T2, z::T3) where {T1, T2, T3}
+function drummond1F1(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
     T = promote_type(T1, T2, T3)
-    absα = T(abs(α))
+    absα = abs(T(α))
     if norm(z) < eps(real(T)) || norm(α) < eps(absα)
         return one(T)
     end
@@ -180,7 +180,7 @@ function drummond1F1(α::T1, β::T2, z::T3) where {T1, T2, T3}
     Nhi /= α+k+1
     Dhi /= α+k+1
     k += 1
-    while abs(Thi-Tmid) > 10*abs(Thi)*eps(real(T)) && abs(Tmid-Tlo) > 10*abs(Tmid)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
         Nhi, Nmid, Nlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Nhi + k*((β+2k+2)*ζ-1)*Nmid + k*(k-1)*ζ*Nlo, Nhi, Nmid
         Dhi, Dmid, Dlo = ((β+k+1)*(k+2)*ζ-(α+2k+1))*Dhi + k*((β+2k+2)*ζ-1)*Dmid + k*(k-1)*ζ*Dlo, Dhi, Dmid
         Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
@@ -195,7 +195,7 @@ function drummond1F1(α::T1, β::T2, z::T3) where {T1, T2, T3}
 end
 
 # ₀F₂(α,β;z)
-function drummond0F2(α::T1, β::T2, z::T3) where {T1, T2, T3}
+function drummond0F2(α::T1, β::T2, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
     T = promote_type(T1, T2, T3)
     if norm(z) < eps(real(T)) || norm(α) < eps(real(T)) || norm(β) < eps(real(T))
         return one(T)
@@ -221,7 +221,7 @@ function drummond0F2(α::T1, β::T2, z::T3) where {T1, T2, T3}
     Dhi, Dmid1, Dmid2, Dlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dmid1 + ζ*k*(k-1)*(3k+α+β+1)*Dmid2 + ζ*k*(k-1)*(k-2)*Dlo, Dhi, Dmid1, Dmid2
     Thi, Tmid1, Tmid2, Tlo = Nhi/Dhi, Thi, Tmid1, Tmid2
     k += 1
-    while abs(Thi-Tmid1) > 10*abs(Thi)*eps(real(T)) && abs(Tmid1-Tmid2) > 10*abs(Tmid1)*eps(real(T)) && abs(Tmid2-Tlo) > 10*abs(Tmid2)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tmid1, Thi, 10eps(real(T)))
         Nhi, Nmid1, Nmid2, Nlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Nhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Nmid1 + ζ*k*(k-1)*(3k+α+β+1)*Nmid2 + ζ*k*(k-1)*(k-2)*Nlo, Nhi, Nmid1, Nmid2
         Dhi, Dmid1, Dmid2, Dlo = (ζ*(k+2)*(α+k+1)*(β+k+1)-1)*Dhi + ζ*k*((k+1)*(α+β+2k)+(α+k)*(β+k)+α+β+3k+2)*Dmid1 + ζ*k*(k-1)*(3k+α+β+1)*Dmid2 + ζ*k*(k-1)*(k-2)*Dlo, Dhi, Dmid1, Dmid2
         Thi, Tmid1, Tmid2, Tlo = Nhi/Dhi, Thi, Tmid1, Tmid2
@@ -231,10 +231,10 @@ function drummond0F2(α::T1, β::T2, z::T3) where {T1, T2, T3}
 end
 
 # ₂F₁(α,β,γ;z)
-function drummond2F1(α::T1, β::T2, γ::T3, z::T4) where {T1, T2, T3, T4}
+function drummond2F1(α::T1, β::T2, γ::T3, z::T4; kmax::Int = 10_000) where {T1, T2, T3, T4}
     T = promote_type(T1, T2, T3, T4)
-    absα = T(abs(α))
-    absβ = T(abs(β))
+    absα = abs(T(α))
+    absβ = abs(T(β))
     if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ)
         return one(T)
     end
@@ -268,7 +268,7 @@ function drummond2F1(α::T1, β::T2, γ::T3, z::T4) where {T1, T2, T3, T4}
     Nhi /= (α+k+1)*(β+k+1)
     Dhi /= (α+k+1)*(β+k+1)
     k += 1
-    while abs(Thi-Tmid) > 10*abs(Thi)*eps(real(T)) && abs(Tmid-Tlo) > 10*abs(Tmid)*eps(real(T)) && k < 10_000
+    while k < kmax && errcheck(Tmid, Thi, 10eps(real(T)))
         Nhi, Nmid, Nlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Nhi + k*((γ+2k+2)*ζ-(α+β+3k))*Nmid + k*(k-1)*(ζ-1)*Nlo, Nhi, Nmid
         Dhi, Dmid, Dlo = ((k+2)*(γ+k+1)*ζ-(α+k+1)*(β+k+1)-k*(α+β+2k+1))*Dhi + k*((γ+2k+2)*ζ-(α+β+3k))*Dmid + k*(k-1)*(ζ-1)*Dlo, Dhi, Dmid
         Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid
@@ -283,9 +283,9 @@ function drummond2F1(α::T1, β::T2, γ::T3, z::T4) where {T1, T2, T3, T4}
 end
 
 # ₘFₙ(α;β;z)
-function drummondpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) where {T1, T2, T3}
+function pFqdrummond(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
     T = promote_type(T1, T2, T3)
-    absα = T.(abs.(α))
+    absα = abs.(T.(α))
     if norm(z) < eps(real(T)) || norm(prod(α)) < eps(prod(absα))
         return one(T)
     end
@@ -295,13 +295,15 @@ function drummondpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) wher
     r = max(p+1, q+2)
     C = zeros(T, r)
     C[1] = one(T)
+    Ĉ = zeros(T, r)
+    Ĉ[1] = one(T)
     P = zeros(T, p+1)
     t = one(T)
     for j in 1:p
         t *= α[j]+1
     end
     P[1] = t
-    err = one(T)
+    err = one(real(T))
     for j in 1:p
         err *= absα[j]+1
     end
@@ -318,7 +320,7 @@ function drummondpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) wher
     D[r+1] = prod(β)*ζ/prod(α)
     R[r+1] = N[r+1]/D[r+1]
     k = 0
-    while (abs(R[r+1]-R[r]) > 10*abs(R[r+1])*eps(real(T)) && k < 10_000) || k < r
+    while k < r || (k < kmax && errcheck(R[r], R[r+1], 10eps(real(T))))
         for j in 1:r
             N[j] = N[j+1]
             D[j] = D[j+1]
@@ -326,25 +328,29 @@ function drummondpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) wher
         end
         t1 = zero(T)
         for j in 0:min(k, q+1)
-            t1 += C[j+1]*Q[j+1]*N[r-j]
+            t1 += Ĉ[j+1]*Q[j+1]*N[r-j]
         end
         if k ≤ q+1
-            t1 += Q[k+1]
+            if p > q
+                t1 += Q[k+1]
+            else
+                t1 += Q[k+1] / T(factorial(k+1))
+            end
         end
         t2 = zero(T)
-        t2 += P[1]*N[r]
+        t2 += Ĉ[1]*P[1]*N[r]
         for j in 1:min(k, p)
-            t2 += C[j+1]*P[j+1]*(N[r-j+1]+N[r-j])
+            t2 += P[j+1]*(C[j+1]*N[r-j+1] + Ĉ[j+1]*N[r-j])
         end
         N[r+1] = ζ*t1-t2
         t1 = zero(T)
         for j in 0:min(k, q+1)
-            t1 += C[j+1]*Q[j+1]*D[r-j]
+            t1 += Ĉ[j+1]*Q[j+1]*D[r-j]
         end
         t2 = zero(T)
-        t2 += P[1]*D[r]
+        t2 += Ĉ[1]*P[1]*D[r]
         for j in 1:min(k, p)
-            t2 += C[j+1]*P[j+1]*(D[r-j+1]+D[r-j])
+            t2 += P[j+1]*(C[j+1]*D[r-j+1] + Ĉ[j+1]*D[r-j])
         end
         D[r+1] = ζ*t1-t2
         R[r+1] = N[r+1]/D[r+1]
@@ -354,14 +360,23 @@ function drummondpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) wher
         N[r+1] /= P[1]
         D[r+1] /= P[1]
         k += 1
-        for j in min(k, max(p, q+1)):-1:1
-            C[j+1] += C[j]
+        if p > q
+            for j in min(k, max(p, q+1)):-1:1
+                C[j+1] += C[j]
+                Ĉ[j+1] += Ĉ[j]
+            end
+        else
+            for j in min(k, max(p, q+1)):-1:1
+                C[j+1] = (C[j+1]*(k+1-j) + C[j])/(k+1)
+                Ĉ[j+1] = (Ĉ[j+1]*(k-j) + Ĉ[j])/(k+1)
+            end
+            Ĉ[1] = Ĉ[1]*k/(k+1)
         end
         t = one(T)
         for j in 1:p
             t *= α[j]+k+1
         end
-        err = one(T)
+        err = one(real(T))
         for j in 1:p
             err *= absα[j]+k+1
         end
diff --git a/src/gauss.jl b/src/gauss.jl
index 9a850d2..007d644 100644
--- a/src/gauss.jl
+++ b/src/gauss.jl
@@ -1,203 +1,145 @@
 # The references to special cases are to Table of Integrals, Series, and Products, § 9.121, followed by NIST's DLMF.
 
-reverseorientation(x::Number) = x
 """
-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)
-  if real(a) > real(b)
-    return _₂F₁(b, a, c, z) # ensure a ≤ b
-  elseif isequal(a, c) # 1. 15.4.6
-    return exp(-b*log1p(-z))
-  elseif isequal(b, c) # 1. 15.4.6
-    return exp(-a*log1p(-z))
-  elseif isequal(c, 0.5)
-    if isequal(a+b, 0) # 31. 15.4.11 & 15.4.12
-      return cosnasinsqrt(2b, z)
-    elseif isequal(a+b, 1) # 32. 15.4.13 & 15.4.14
-      return cosnasinsqrt(1-2b, z)*exp(-0.5log1p(-z))
-    elseif isequal(b-a, 0.5) # 15.4.7 & 15.4.8
-      return expnlog1pcoshatanhsqrt(-2a, z)
+    if real(b) < real(a)
+        return _₂F₁(b, a, c, z) # ensure a ≤ b
+    elseif isequal(a, c) # 1. 15.4.6
+        return exp(-b*log1p(-z))
+    elseif isequal(b, c) # 1. 15.4.6
+        return exp(-a*log1p(-z))
+    elseif isequal(c, 0.5)
+        if isequal(a+b, 0) # 31. 15.4.11 & 15.4.12
+            return cosnasinsqrt(2b, z)
+        elseif isequal(a+b, 1) # 32. 15.4.13 & 15.4.14
+            return cosnasinsqrt(1-2b, z)*exp(-0.5log1p(-z))
+        elseif isequal(b-a, 0.5) # 15.4.7 & 15.4.8
+            return expnlog1pcoshatanhsqrt(-2a, z)
+        end
+    elseif isequal(c, 1.5)
+        if abeqcd(a, b, 0.5) # 13. 15.4.4 & 15.4.5
+            return sqrtasinsqrt(z)
+        elseif abeqcd(a, b, 1) # 14.
+            return sqrtasinsqrt(z)*exp(-0.5log1p(-z))
+        elseif abeqcd(a, b, 0.5, 1) # 15. 15.4.2 & 15.4.3
+            return sqrtatanhsqrt(z)
+        elseif isequal(a+b, 1) # 29. 15.4.15 & 15.4.16
+            return sinnasinsqrt(1-2b, z)
+        elseif isequal(a+b, 2) # 30.
+            return sinnasinsqrt(2-2b, z)*exp(-0.5log1p(-z))
+        elseif isequal(b-a, 0.5) # 4. 15.4.9 & 15.4.10
+            return expnlog1psinhatanhsqrt(1-2a, z)
+        end
+    elseif isequal(c, 2)
+        if abeqcd(a, b, 1) # 6. 15.4.1
+            return log1pover(-z)
+        elseif a ∈ ℤ && b == 1 # 5.
+            return expm1nlog1p(1-a, -z)
+        elseif a == 1 && b ∈ ℤ # 5.
+            return expm1nlog1p(1-b, -z)
+        end
+    elseif isequal(c, 4)
+        if abeqcd(a, b, 2)
+            return logandpoly(z)
+        end
+    elseif isequal(c, 2.5) && abeqcd(a, b, 1, 1.5)
+        return speciallog(z)
     end
-  elseif isequal(c, 1.5)
-    if abeqcd(a, b, 0.5) # 13. 15.4.4 & 15.4.5
-      return sqrtasinsqrt(z)
-    elseif abeqcd(a, b, 1) # 14.
-      return sqrtasinsqrt(z)*exp(-0.5log1p(-z))
-    elseif abeqcd(a, b, 0.5, 1) # 15. 15.4.2 & 15.4.3
-      return sqrtatanhsqrt(z)
-    elseif isequal(a+b, 1) # 29. 15.4.15 & 15.4.16
-      return sinnasinsqrt(1-2b, z)
-    elseif isequal(a+b, 2) # 30.
-      return sinnasinsqrt(2-2b, z)*exp(-0.5log1p(-z))
-    elseif isequal(b-a, 0.5) # 4. 15.4.9 & 15.4.10
-      return expnlog1psinhatanhsqrt(1-2a, z)
-    end
-  elseif isequal(c, 2)
-    if abeqcd(a, b, 1) # 6. 15.4.1
-      return log1pover(-z)
-    elseif a ∈ ℤ && b == 1 # 5.
-      return expm1nlog1p(1-a, -z)
-    elseif a == 1 && b ∈ ℤ # 5.
-      return expm1nlog1p(1-b, -z)
-    end
-  elseif isequal(c, 4)
-    if abeqcd(a, b, 2)
-      return logandpoly(z)
-    end
-  elseif isequal(c, 2.5) && abeqcd(a, b, 1, 1.5)
-    return speciallog(z)
-  end
-  _₂F₁general(a, b, c, z) # catch-all
+    return _₂F₁general(a, b, c, z) # catch-all
 end
 
-
-
-
 # Special case of (-x)^a*_₂F₁ to handle LogNumber correctly in RiemannHilbert.jl
 function mxa_₂F₁(a, b, c, z)
-  if isequal(c, 2)
-    if abeqcd(a, b, 1) # 6. 15.4.1
-      return log1p(-z)
-    end
-  elseif isequal(c, 4)
-    if abeqcd(a, b, 2)
-      return 6*(-2 + (1-2/z)*log1p(-z))
+    if isequal(c, 2)
+        if abeqcd(a, b, 1) # 6. 15.4.1
+            return log1p(-z)
+        end
+    elseif isequal(c, 4)
+        if abeqcd(a, b, 2)
+            return 6*(-2 + (1-2/z)*log1p(-z))
+        end
     end
-  end
-  (-z)^a*_₂F₁(a, b, c, z)
+    return (-z)^a*_₂F₁(a, b, c, z)
 end
 
-
-_₂F₁(a::Number, b::Number, c::Number, z::AbstractArray) = reshape(promote_type(typeof(a), typeof(b), typeof(c), eltype(z))[ _₂F₁(a, b, c, z[i]) for i in eachindex(z) ], size(z))
-
 """
-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.
 """
 function _₂F₁general(a, b, c, z)
-  T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
-
-  real(b) < real(a) && return _₂F₁general(b, a, c, z)
-  real(c) < real(a) + real(b) && return exp((c-a-b)*log1p(-z))*_₂F₁general(c-a, c-b, c, z)
-
-  if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
-    _₂F₁maclaurin(a, b, c, z)
-  elseif abs(z/(z-1)) ≤ ρ
-    exp(-a*log1p(-z))_₂F₁maclaurin(a, c-b, c, z/(z-1))
-  elseif abs(inv(z)) ≤ ρ
-    _₂F₁Inf(a, b, c, z)
-  elseif abs(1-inv(z)) ≤ ρ
-    exp(-a*log1p(-z))*_₂F₁Inf(a, c-b, c, reverseorientation(z/(z-1)))
-  elseif abs(1-z) ≤ ρ
-    _₂F₁one(a, b, c, z)
-  elseif abs(inv(1-z)) ≤ ρ
-    exp(-a*log1p(-z))*_₂F₁one(a, c-b, c, reverseorientation(z/(z-1)))
-  else
-    _₂F₁taylor(a, b, c, z)
-  end
+    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
+
+    real(b) < real(a) && return _₂F₁general(b, a, c, z)
+    real(c) < real(a) + real(b) && return exp((c-a-b)*log1p(-z))*_₂F₁general(c-a, c-b, c, z)
+
+    if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
+        return _₂F₁maclaurin(a, b, c, z)
+    elseif abs(z/(z-1)) ≤ ρ
+        return exp(-a*log1p(-z))_₂F₁maclaurin(a, c-b, c, z/(z-1))
+    elseif abs(inv(z)) ≤ ρ
+        return _₂F₁Inf(a, b, c, z)
+    elseif abs(1-inv(z)) ≤ ρ
+        return exp(-a*log1p(-z))*_₂F₁Inf(a, c-b, c, z/(z-1))
+    elseif abs(1-z) ≤ ρ
+        return _₂F₁one(a, b, c, z)
+    elseif abs(inv(1-z)) ≤ ρ
+        return exp(-a*log1p(-z))*_₂F₁one(a, c-b, c, z/(z-1))
+    else
+        return pFqweniger([a, b], [c], z) # _₂F₁taylor(a, b, c, z)
+    end
 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.
 """
 function _₂F₁general2(a, b, c, z)
-  T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
-  if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
-    _₂F₁maclaurin(a, b, c, z)
-  elseif abs(z / (z - 1)) ≤ ρ && absarg(1 - z) < convert(real(T), π) # 15.8.1
-    w = z/(z-1)
-    _₂F₁maclaurin(a, c-b, c, w)*exp(-a*log1p(-z))
-  elseif abs(inv(z)) ≤ ρ && absarg(-z) < convert(real(T), π)
-    w = inv(z)
-    if isapprox(a, b) # 15.8.8
-      gamma(c)/gamma(a)/gamma(c-a)*(-w)^a*_₂F₁logsumalt(a, c-a, z, w)
-    elseif a-b ∉ ℤ # 15.8.2
-      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))
-    else
-      # TODO: full 15.8.8
-      wenigerpFq([a, b], [c], z)
-    end
-  elseif abs(inv(1-z)) ≤ ρ && absarg(-z) < convert(real(T), π)
-    w = inv(1-z)
-    if isapprox(a, b) # 15.8.9
-      gamma(c)*exp(-a*log1p(-z))/gamma(a)/gamma(c-b)*_₂F₁logsum(a, c-a, z, w, 1)
-    elseif a-b ∉ ℤ # 15.8.3
-      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))
-    else
-      # TODO: full 15.8.9
-      wenigerpFq([a, b], [c], z)
-    end
-  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
-      gamma(c)/gamma(a)/gamma(b)*_₂F₁logsum(a, b, z, w, -1)
-    elseif c - a - b ∉ ℤ # 15.8.4
-      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))
-    else
-      # TODO: full 15.8.10
-      wenigerpFq([a, b], [c], z)
+    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
+    if abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
+        return _₂F₁maclaurin(a, b, c, z)
+    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))
+    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)
+        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))
+        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)
+        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))
+        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)
+        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))
+        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)
+        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))
+        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)
+        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))
+        end
     end
-  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
-      gamma(c)*z^(-a)/gamma(a)*_₂F₁logsumalt(a, b, z, w)
-    elseif c - a - b ∉ ℤ # 15.8.5
-      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))
-    else
-      # TODO: full 15.8.11
-      wenigerpFq([a, b], [c], z)
-    end
-  elseif abs(z-0.5) > 0.5
-    if isapprox(a, b) && !isapprox(c, a+0.5)
-      gamma(c)/gamma(a)/gamma(c-a)*(0.5-z)^(-a)*_₂F₁continuationalt(a, c, 0.5, z)
-    elseif a-b ∉ ℤ
-      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))
-    else
-      wenigerpFq([a, b], [c], z)
-    end
-  else
-    throw(DomainError())
-  end
-end
-
-"""
-Compute the generalized hypergeometric function `₃F₂(a₁, a₂, a₃;b₁, b₂;z)` with the Maclaurin series.
-"""
-function _₃F₂(a₁, a₂, a₃, b₁, b₂, z)
-  if abs(z) ≤ ρ
-    _₃F₂maclaurin(a₁, a₂, a₃, b₁, b₂, z)
-  else
-    wenigerpFq([a₁, a₂, a₃], [b₁, b₂], z)
-  end
-end
-_₃F₂(a₁, b₁, z) = _₃F₂(1, 1, a₁, 2, b₁, z)
-
-"""
-Compute the generalized hypergeometric function `mFn(a;b;z)` with the Maclaurin series.
-"""
-function mFn(a::AbstractVector, b::AbstractVector, z)
-    if abs(z) ≤ ρ || length(a) ≤ length(b)
-        mFnmaclaurin(a, b, z)
-    else
-        wenigerpFq(a, b, z)
-    end
-end
-
-"""
-computation of the generalized hypergeometric function `mFn(a;b;z) by continued
-fraction`
-"""
-function mFncontinuedfraction(a::AbstractVector{S}, b::AbstractVector{U}, z::V) where {S, U, V}
-  T = promote_type(S, U, V)
-  numerator(i) = - z * prod(i .+ a) / prod(i .+ b) / (i + 1)
-  denominator(i) = 1 - numerator(i)
-  K = continuedfraction(denominator, numerator, 10eps(real(T))) - denominator(0)
-  # iszero(K + 1) leads to Inf, when e.g. Float64s run out of digits
-  return 1 + z * prod(a) / prod(b) / (1 + K)
+    return pFqweniger([a, b], [c], z)
 end
diff --git a/src/generalized.jl b/src/generalized.jl
new file mode 100644
index 0000000..bdacc3c
--- /dev/null
+++ b/src/generalized.jl
@@ -0,0 +1,44 @@
+@doc raw"""
+    pFq(α, β; z)
+Compute the generalized hypergeometric function, defined by
+```math
+\operatorname{pFq}(α, β; z) = \sum_{k=0}^\infty \dfrac{(\alpha_1)_k\cdots(\alpha_p)_k}{(\beta_1)_k\cdots(\beta_q)_k}\dfrac{z^k}{k!}.
+```
+"""
+function pFq(α::AbstractVector, β::AbstractVector, z; kwds...)
+    if length(α) == length(β) == 0
+        return exp(z)
+    elseif length(α) == 1 && length(β) == 0
+        return exp(-α[1]*log1p(-z))
+    elseif length(α) == 2 && length(β) == 1
+        return _₂F₁(α[1], α[2], β[1], float(z))
+    elseif abs(z) ≤ ρ || length(α) ≤ length(β)
+        return pFqmaclaurin(α, β, float(z); kwds...)
+    else
+        return pFqweniger(α, β, float(z); kwds...)
+    end
+end
+
+"""
+Compute the generalized hypergeometric function `pFq(α, β; z)` by continued fraction.
+"""
+function pFqcontinuedfraction(α::AbstractVector{S}, β::AbstractVector{U}, z::V) where {S, U, V}
+    T = promote_type(S, U, V)
+    numerator(i) = - z * prod(i .+ α) / prod(i .+ β) / (i + 1)
+    denominator(i) = 1 - numerator(i)
+    K = continuedfraction(denominator, numerator, 10eps(real(T))) - denominator(0)
+    # iszero(K + 1) leads to Inf, when e.g. Float64s run out of digits
+    return 1 + z * prod(α) / prod(β) / (1 + K)
+end
+
+"""
+Compute the generalized hypergeometric function `₃F₂(a₁, a₂, a₃, b₁, b₂; z)`.
+"""
+function _₃F₂(a₁, a₂, a₃, b₁, b₂, z; kwds...)
+    if abs(z) ≤ ρ
+        _₃F₂maclaurin(a₁, a₂, a₃, b₁, b₂, float(z); kwds...)
+    else
+        pFqweniger([a₁, a₂, a₃], [b₁, b₂], float(z); kwds...)
+    end
+end
+_₃F₂(a₁, b₁, z; kwds...) = _₃F₂(1, 1, a₁, 2, b₁, z; kwds...)
diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl
index 10c748c..4b4963f 100644
--- a/src/specialfunctions.jl
+++ b/src/specialfunctions.jl
@@ -1,5 +1,4 @@
-import DualNumbers: Dual, Dual128, DualComplex256
-import SpecialFunctions: gamma, digamma, factorial
+@inline errcheck(x, y, tol) = isfinite(x) && isfinite(y) && (norm(x-y) > max(norm(x), norm(y))*tol)
 
 # Same as in FastTransforms.jl
 """
@@ -31,17 +30,22 @@ end
 
 ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))
 
+"""
+    @clenshaw(x, c...)
+
+Compute `∑_{i=1}^N cᵢ Tᵢ₋₁(x)`.
+"""
 macro clenshaw(x, c...)
     a, b = :(zero(t)), :(zero(t))
     as = []
     N = length(c)
     for k = N:-1:2
-        ak = Symbol("a",k)
+        ak = Symbol("a", k)
         push!(as, :($ak = $a))
-        a = :(muladd(t,$a,$(esc(c[k]))-$b))
+        a = :(muladd(t, $a, $(esc(c[k]))-$b))
         b = :($ak)
     end
-    ex = Expr(:block,as...,:(muladd(t/2,$a,$(esc(c[1]))-$b)))
+    ex = Expr(:block, as..., :(muladd(t/2, $a, $(esc(c[1]))-$b)))
     Expr(:block, :(t = $(esc(2))*$(esc(x))), ex)
 end
 
@@ -94,11 +98,11 @@ expm1nlog1p(n, x) = x == 0 ? one(x) : expm1(n*log1p(x))/(n*x)
 log1pover(s) = log1p(s)/s
 logandpoly(x) = x == 0 ? one(x) : 6*(-2x+(x-2)*log1p(-x))/x^3
 function logandpoly(x::Union{Float64, ComplexF64})
-  if abs(x) > 0.2
-    6*(-2x+(x-2)*log1p(-x))/x^3
-  else
-    logandpolyseries(x)
-  end
+    if abs(x) > 0.2
+        6*(-2x+(x-2)*log1p(-x))/x^3
+    else
+        logandpolyseries(x)
+    end
 end
 
 logandpolyseries(x::Union{Float64, Dual128, ComplexF64, DualComplex256}) = @evalpoly(x, 1.0, 1.0, 0.9, 0.8, 0.7142857142857143, 0.6428571428571429, 0.5833333333333334, 0.5333333333333333, 0.4909090909090909, 0.45454545454545453, 0.4230769230769231, 0.3956043956043956, 0.37142857142857144, 0.35, 0.33088235294117646, 0.3137254901960784, 0.2982456140350877, 0.28421052631578947, 0.2714285714285714, 0.2597402597402597)
@@ -106,23 +110,23 @@ logandpolyseries(x::Union{Float64, Dual128, ComplexF64, DualComplex256}) = @eval
 speciallog(x::Real) = iszero(x) ? one(x) : (x > 0 ? (s = sqrt(x); 3(atanh(s)-s)/s^3) : (s = sqrt(-x); 3(s-atan(s))/s^3))
 speciallog(x) = iszero(x) ? one(x) : (s = sqrt(-x); 3(s-atan(s))/s^3)
 function speciallog(x::Float64)
-  if x > 0.2
-    s = sqrt(x)
-    3(atanh(s)-s)/s^3
-  elseif x < -0.2
-    s = sqrt(-x)
-    3(s-atan(s))/s^3
-  else
-    speciallogseries(x)
-  end
+    if x > 0.2
+        s = sqrt(x)
+        3(atanh(s)-s)/s^3
+    elseif x < -0.2
+        s = sqrt(-x)
+        3(s-atan(s))/s^3
+    else
+        speciallogseries(x)
+    end
 end
 function speciallog(x::ComplexF64)
-  if abs(x) > 0.2
-    s = sqrt(-x)
-    3(s-atan(s))/s^3
-  else
-    speciallogseries(x)
-  end
+    if abs(x) > 0.2
+        s = sqrt(-x)
+        3(s-atan(s))/s^3
+    else
+        speciallogseries(x)
+    end
 end
 # The Maclaurin series fails to be accurate to 1e-15 near x ≈ ±0.2. So we use a highly accurate Chebyshev expansion.
 speciallogseries(x::Union{Float64, Dual128}) = @clenshaw(5.0x, 1.0087391788544393911192, 1.220474262857857637288e-01, 8.7957928919918696061703e-03, 6.9050958578444820505037e-04, 5.7037120050065804396306e-05, 4.8731405131379353370205e-06, 4.2648797509486828820613e-07, 3.800372208946157617901e-08, 3.434168059359993493634e-09, 3.1381484326392473547608e-10, 2.8939845618385022798906e-11, 2.6892186934806386106143e-12, 2.5150879096374730760324e-13, 2.3652490233687788117887e-14, 2.2349973917002118259929e-15, 2.120769988408948118084e-16)
@@ -136,83 +140,83 @@ unsafe_gamma(x::Float64) = ccall((:tgamma, libm),  Float64, (Float64, ), x)
 unsafe_gamma(x::Float32) = ccall((:tgammaf, libm),  Float32, (Float32, ), x)
 unsafe_gamma(x::Real) = unsafe_gamma(float(x))
 function unsafe_gamma(x::BigFloat)
-  z = BigFloat()
-  ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
-  return z
+    z = BigFloat()
+    ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
+    return z
 end
 unsafe_gamma(z::Complex) = gamma(z)
 unsafe_gamma(z::Dual) = (r = realpart(z);w = unsafe_gamma(r); dual(w, w*digamma(r)*dualpart(z)))
 
 """
-`@lanczosratio(z, ϵ, c₀, c...)`
+    @lanczosratio(z, ϵ, c₀, c...)
 
-Compute `∑_{i=1}^N cᵢ/(z-1+i)/(z-1+i+ϵ) / ( c₀ + ∑_{i=1}^N cᵢ/(z-1+i) )`
+Compute `∑_{i=1}^N cᵢ/(z-1+i)/(z-1+i+ϵ) / ( c₀ + ∑_{i=1}^N cᵢ/(z-1+i) )`.
 """
 macro lanczosratio(z, ϵ, c₀, c...)
-  ex_num = :(zero(zm1))
-  ex_den = esc(c₀)
-  for i = 0:length(c)-1
-    temp = :(inv(zm1+$i))
-    ex_num = :(muladd($(esc(c[i+1])), $temp*inv(zm1pϵ+$i), $ex_num))
-    ex_den = :(muladd($(esc(c[i+1])), $temp, $ex_den))
-  end
-  ex = :($ex_num/$ex_den)
-  Expr(:block, :(zm1 = $(esc(z))), :(zm1pϵ = $(esc(z))+$(esc(ϵ))), ex)
+    ex_num = :(zero(zm1))
+    ex_den = esc(c₀)
+    for i = 0:length(c)-1
+        temp = :(inv(zm1+$i))
+        ex_num = :(muladd($(esc(c[i+1])), $temp*inv(zm1pϵ+$i), $ex_num))
+        ex_den = :(muladd($(esc(c[i+1])), $temp, $ex_den))
+    end
+    ex = :($ex_num/$ex_den)
+    Expr(:block, :(zm1 = $(esc(z))), :(zm1pϵ = $(esc(z))+$(esc(ϵ))), ex)
 end
 
 lanczosratio(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256}) = @lanczosratio(z, ϵ, 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699E-4, 0.46523628927048575665E-4, -0.98374475304879564677E-4, 0.15808870322491248884E-3, -0.21026444172410488319E-3, 0.21743961811521264320E-3, -0.16431810653676389022E-3, 0.84418223983852743293E-4, -0.26190838401581408670E-4, 0.36899182659531622704E-5)
 
 function lanczosapprox(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
-  zm0p5 = z-0.5
-  zpgm0p5 = zm0p5+4.7421875
-  return zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ))
+    zm0p5 = z-0.5
+    zpgm0p5 = zm0p5+4.7421875
+    return zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ))
 end
 
 function H(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
-  zm0p5 = z-0.5
-  zpgm0p5 = zm0p5+4.7421875
-  if real(z) ≥ 1/2
-    if z == z+ϵ # ϵ is numerical 0
-      zm0p5/zpgm0p5 + log(zpgm0p5) - 1 - lanczosratio(z, ϵ)
-    else
-      expm1( zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ)) )/ϵ
-    end
-  else
-    tpz = tanpi(z)
-    if z == z+ϵ # ϵ is numerical 0
-      H(1-z, ϵ) - π/tpz
+    zm0p5 = z-0.5
+    zpgm0p5 = zm0p5+4.7421875
+    if real(z) ≥ 1/2
+        if z == z+ϵ # ϵ is numerical 0
+            zm0p5/zpgm0p5 + log(zpgm0p5) - 1 - lanczosratio(z, ϵ)
+        else
+            expm1( zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ)) )/ϵ
+        end
     else
-      temp = (cospi(ϵ) + sinpi(ϵ)/tpz)*H(1-z, -ϵ) + .5ϵ*(π*sinc(.5ϵ))^2 - π*sinc(ϵ)/tpz
-      temp/(1-ϵ*temp)
+        tpz = tanpi(z)
+        if z == z+ϵ # ϵ is numerical 0
+            H(1-z, ϵ) - π/tpz
+        else
+            temp = (cospi(ϵ) + sinpi(ϵ)/tpz)*H(1-z, -ϵ) + .5ϵ*(π*sinc(.5ϵ))^2 - π*sinc(ϵ)/tpz
+            temp/(1-ϵ*temp)
+        end
     end
-  end
 end
 
 """
 Compute the function (1/Γ(z)-1/Γ(z+ϵ))/ϵ
 """
 function G(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
-  n, zpϵ = round(Int, real(z)), z+ϵ
-  if abs(ϵ) > 0.1
-    (inv(unsafe_gamma(z))-inv(unsafe_gamma(zpϵ)))/ϵ
-  elseif z ≠ zpϵ
-    m = round(Int, real(zpϵ))
-    if z == n && n ≤ 0
-      -inv(ϵ*unsafe_gamma(zpϵ))
-    elseif zpϵ == m && m ≤ 0
-      inv(ϵ*unsafe_gamma(z))
-    elseif abs(z+abs(n)) < abs(zpϵ+abs(m))
-      H(z, ϵ)/unsafe_gamma(zpϵ)
-    else
-      H(zpϵ, -ϵ)/unsafe_gamma(z)
-    end
-  else # ϵ is numerical 0
-    if z == n && n ≤ 0
-      (-1)^(n+1)*unsafe_gamma(1-n)
-    else
-      digamma(z)/unsafe_gamma(z)
+    n, zpϵ = round(Int, real(z)), z+ϵ
+    if abs(ϵ) > 0.1
+        (inv(unsafe_gamma(z))-inv(unsafe_gamma(zpϵ)))/ϵ
+    elseif z ≠ zpϵ
+        m = round(Int, real(zpϵ))
+        if z == n && n ≤ 0
+            -inv(ϵ*unsafe_gamma(zpϵ))
+        elseif zpϵ == m && m ≤ 0
+            inv(ϵ*unsafe_gamma(z))
+        elseif abs(z+abs(n)) < abs(zpϵ+abs(m))
+            H(z, ϵ)/unsafe_gamma(zpϵ)
+        else
+            H(zpϵ, -ϵ)/unsafe_gamma(z)
+        end
+    else # ϵ is numerical 0
+        if z == n && n ≤ 0
+            (-1)^(n+1)*unsafe_gamma(1-n)
+        else
+            digamma(z)/unsafe_gamma(z)
+        end
     end
-  end
 end
 
 G(z::Number, ϵ::Number) = ϵ == 0 ? digamma(z)/unsafe_gamma(z) : (inv(unsafe_gamma(z))-inv(unsafe_gamma(z+ϵ)))/ϵ
@@ -222,42 +226,42 @@ G(z::Number, ϵ::Number) = ϵ == 0 ? digamma(z)/unsafe_gamma(z) : (inv(unsafe_ga
 Compute the function ((z+ϵ)ₘ-(z)ₘ)/ϵ
 """
 function P(z::Number, ϵ::Number, m::Int)
-  n₀ = -round(Int, real(z))
-  if ϵ == 0
-    if 0 ≤ n₀ < m
-      ret1, ret2, n = one(z), zero(z), 0
-      while n < m
-        n == n₀ && (n+=1; continue)
-        ret1 *= z+n
-        ret2 += inv(z+n)
-        n+=1
-      end
-      ret1 + pochhammer(z, m)*ret2
-    else
-      ret = zero(z)
-      for n=0:m-1
-        ret += inv(z+n)
-      end
-      pochhammer(z, m)*ret
-    end
-  else
-    if 0 ≤ n₀ < m
-      zpϵ, ret1, ret2, n = z+ϵ, one(z), zero(z), 0
-      while n < m
-        n == n₀ && (n+=1; continue)
-        ret1 *= zpϵ+n
-        ret2 += log1p(ϵ/(z+n))
-        n+=1
-      end
-      ret1 + pochhammer(z, m)*expm1(ret2)/ϵ
+    n₀ = -round(Int, real(z))
+    if ϵ == 0
+        if 0 ≤ n₀ < m
+            ret1, ret2, n = one(z), zero(z), 0
+            while n < m
+                n == n₀ && (n+=1; continue)
+                ret1 *= z+n
+                ret2 += inv(z+n)
+                n += 1
+            end
+            ret1 + pochhammer(z, m)*ret2
+        else
+            ret = zero(z)
+            for n=0:m-1
+                ret += inv(z+n)
+            end
+            pochhammer(z, m)*ret
+        end
     else
-      ret = zero(z)
-      for n=0:m-1
-        ret += log1p(ϵ/(z+n))
-      end
-      pochhammer(z, m)*expm1(ret)/ϵ
+        if 0 ≤ n₀ < m
+            zpϵ, ret1, ret2, n = z+ϵ, one(z), zero(z), 0
+            while n < m
+                n == n₀ && (n+=1; continue)
+                ret1 *= zpϵ+n
+                ret2 += log1p(ϵ/(z+n))
+                n += 1
+            end
+            ret1 + pochhammer(z, m)*expm1(ret2)/ϵ
+        else
+            ret = zero(z)
+            for n=0:m-1
+                ret += log1p(ϵ/(z+n))
+            end
+            pochhammer(z, m)*expm1(ret)/ϵ
+        end
     end
-  end
 end
 
 E(z::Number, ϵ::Number) = ϵ == 0 ? z : expm1(ϵ*z)/ϵ
@@ -271,33 +275,32 @@ reconeβ₀(a, b, c, w, m::Int, ϵ) = abs(ϵ) > 0.1 ? ( pochhammer(float(a), m)*
 reconeγ₀(a, b, c, w, m::Int, ϵ) = gamma(c)*pochhammer(float(a), m)*pochhammer(b, m)*w^m/(gamma(a+m+ϵ)*gamma(b+m+ϵ)*gamma(m+1)*gamma(1-ϵ))
 
 function Aone(a, b, c, w, m::Int, ϵ)
-  αₙ = reconeα₀(a, b, c, m, ϵ)*one(w)
-  ret = m ≤ 0 ? zero(w) : αₙ
-  for n = 0:m-2
-    αₙ *= (a+n)*(b+n)/((n+1)*(1-m-ϵ+n))*w
-    ret += αₙ
-  end
-  ret
+    αₙ = reconeα₀(a, b, c, m, ϵ)*one(w)
+    ret = m ≤ 0 ? zero(w) : αₙ
+    for n = 0:m-2
+        αₙ *= (a+n)*(b+n)/((n+1)*(1-m-ϵ+n))*w
+        ret += αₙ
+    end
+    ret
 end
 
 function Bone(a, b, c, w, m::Int, ϵ)
-  βₙ, γₙ = reconeβ₀(a, b, c, w, m, ϵ)*one(w), reconeγ₀(a, b, c, w, m, ϵ)*w
-  ret, err, n = βₙ, 1.0, 0
-  while err > 10eps()
-    βₙ = (a+m+n+ϵ)*(b+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(b+m+n)/(m+n+1) - (a+m+n) - (b+m+n) - ϵ + (a+m+n+ϵ)*(b+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
-    ret += βₙ
-    err = errcheck(βₙ, ret)
-    γₙ *= (a+m+n)*(b+m+n)/((m+n+1)*(n+1-ϵ))*w
-    n+=1
-  end
-  ret
+    βₙ, γₙ = reconeβ₀(a, b, c, w, m, ϵ)*one(w), reconeγ₀(a, b, c, w, m, ϵ)*w
+    ret, n = βₙ, 0
+    while abs(βₙ) > 10abs(ret)*eps() || n ≤ 0
+        βₙ = (a+m+n+ϵ)*(b+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(b+m+n)/(m+n+1) - (a+m+n) - (b+m+n) - ϵ + (a+m+n+ϵ)*(b+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
+        ret += βₙ
+        γₙ *= (a+m+n)*(b+m+n)/((m+n+1)*(n+1-ϵ))*w
+        n += 1
+    end
+    ret
 end
 
 function _₂F₁one(a, b, c, z)
-  m = round(Int, real(c-(a+b)))
-  ϵ = c-(a+b)-m
-  w = 1-z
-  (-1)^m/sinc(ϵ)*(Aone(a, b, c, w, m, ϵ) + Bone(a, b, c, w, m, ϵ))
+    m = round(Int, real(c-(a+b)))
+    ϵ = c-(a+b)-m
+    w = 1-z
+    (-1)^m/sinc(ϵ)*(Aone(a, b, c, w, m, ϵ) + Bone(a, b, c, w, m, ϵ))
 end
 
 # Transformation formula w = 1/z
@@ -312,226 +315,192 @@ recInfβ₀(a, b, c, w, m::Int, ϵ) = abs(ϵ) > 0.1 ?
 recInfγ₀(a, b, c, w, m::Int, ϵ) = gamma(c)*pochhammer(float(a), m)*pochhammer(float(1-c+a), m)*w^m/(gamma(a+m+ϵ)*gamma(c-a)*gamma(m+1)*gamma(1-ϵ))
 
 function AInf(a, b, c, w, m::Int, ϵ)
-  αₙ = recInfα₀(a, b, c, m, ϵ)*one(w)
-  ret = m ≤ 0 ? zero(w) : αₙ
-  for n = 0:m-2
-    αₙ *= (a+n)*(1-c+a+n)/((n+1)*(1-m-ϵ+n))*w
-    ret += αₙ
-  end
-  ret
+    αₙ = recInfα₀(a, b, c, m, ϵ)*one(w)
+    ret = m ≤ 0 ? zero(w) : αₙ
+    for n = 0:m-2
+        αₙ *= (a+n)*(1-c+a+n)/((n+1)*(1-m-ϵ+n))*w
+        ret += αₙ
+    end
+    ret
 end
 
 function BInf(a, b, c, win, m::Int, ϵ)
-  w=win
-  βₙ, γₙ = recInfβ₀(a, b, c, win, m, ϵ)*one(w), recInfγ₀(a, b, c, win, m, ϵ)*w
-  ret, err, n = βₙ, 1.0, 0
-  while err > 10eps()
-    βₙ = (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(1-c+a+m+n)/(m+n+1) - (a+m+n) - (1-c+a+m+n) - ϵ + (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
-    ret += βₙ
-    err = errcheck(βₙ, ret)
-    γₙ *= (a+m+n)*(1-c+a+m+n)/((m+n+1)*(n+1-ϵ))*w
-    n+=1
-  end
-  ret
+    w = win
+    βₙ, γₙ = recInfβ₀(a, b, c, win, m, ϵ)*one(w), recInfγ₀(a, b, c, win, m, ϵ)*w
+    ret, n = βₙ, 0
+    while abs(βₙ) > 10abs(ret)*eps() || n ≤ 0
+        βₙ = (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(1-c+a+m+n)/(m+n+1) - (a+m+n) - (1-c+a+m+n) - ϵ + (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
+        ret += βₙ
+        γₙ *= (a+m+n)*(1-c+a+m+n)/((m+n+1)*(n+1-ϵ))*w
+        n += 1
+    end
+    ret
 end
 
 function _₂F₁Inf(a, b, c, z)
-  m = round(Int, real(b-a))
-  ϵ = b-a-m
-  w = reverseorientation(inv(z))  # we've swapped the branch cut
-  (-1)^m*(-w)^a/sinc(ϵ)*(AInf(a, b, c, w, m, ϵ) + BInf(a, b, c, w, m, ϵ))
+    m = round(Int, real(b-a))
+    ϵ = b-a-m
+    w = inv(z)
+    (-1)^m*(-w)^a/sinc(ϵ)*(AInf(a, b, c, w, m, ϵ) + BInf(a, b, c, w, m, ϵ))
 end
 
 
 function _₂F₁maclaurin(a::Number, b::Number, c::Number, z::Number)
-  T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
-  S₀, S₁, err, j = one(T), one(T)+a*b*z/c, one(real(T)), 1
-  while err > 10eps(real(T))
-    rⱼ = (a+j)/(j+1)*(b+j)/(c+j)
-    S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
-    err = errcheck(S₁-S₀, S₀)
-    j+=1
-  end
-  return S₁
+    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
+    S₀, S₁, j = one(T), one(T)+a*b*z/c, 1
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
+        rⱼ = (a+j)/(j+1)*(b+j)/(c+j)
+        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
+        j += 1
+    end
+    return S₁
 end
 
 function _₂F₁maclaurinalt(a::Number, b::Number, c::Number, z::Number)
-  T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
-  C, S, err, j = one(T), one(T), one(real(T)), 0
-  while err > 10eps(real(T))
-    C *= (a+j)/(j+1)*(b+j)/(c+j)*z
-    S += C
-    err = errcheck(C, S)
-    j+=1
-  end
-  return S
+    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
+    C, S, j = one(T), one(T), 0
+    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
+        C *= (a+j)/(j+1)*(b+j)/(c+j)*z
+        S += C
+        j += 1
+    end
+    return S
 end
 
 function _₂F₁continuation(s::Number, t::Number, c::Number, z₀::Number, z::Number)
-  T = promote_type(typeof(s), typeof(t), typeof(c), typeof(z₀), typeof(z))
-  izz₀, d0, d1 = inv(z-z₀), one(T), s/(2s-t+one(T))*((s+1)*(1-2z₀)+(t+1)*z₀-c)
-  S₀, S₁, izz₀j, err, j = one(T), one(T)+d1*izz₀, izz₀, one(real(T)), 2
-  while err > 10eps(real(T))
-    d0, d1, izz₀j = d1, (j+s-one(T))/j/(j+2s-t)*(((j+s)*(1-2z₀)+(t+1)*z₀-c)*d1 + z₀*(1-z₀)*(j+s-2)*d0), izz₀j*izz₀
-    S₀, S₁ = S₁, S₁+d1*izz₀j
-    err = errcheck(S₁-S₀, S₀)
-    j+=1
-  end
-  return S₁
+    T = promote_type(typeof(s), typeof(t), typeof(c), typeof(z₀), typeof(z))
+    izz₀, d0, d1 = inv(z-z₀), one(T), s/(2s-t+one(T))*((s+1)*(1-2z₀)+(t+1)*z₀-c)
+    S₀, S₁, izz₀j, j = one(T), one(T)+d1*izz₀, izz₀, 2
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 2
+        d0, d1, izz₀j = d1, (j+s-one(T))/j/(j+2s-t)*(((j+s)*(1-2z₀)+(t+1)*z₀-c)*d1 + z₀*(1-z₀)*(j+s-2)*d0), izz₀j*izz₀
+        S₀, S₁ = S₁, S₁+d1*izz₀j
+        j += 1
+    end
+    return S₁
 end
 
 function _₂F₁continuationalt(a::Number, c::Number, z₀::Number, z::Number)
-  T = promote_type(typeof(a), typeof(c), typeof(z₀), typeof(z))
-  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)
-  S₀ = cⱼ
-  cⱼ += 2/one(T)-one(T)/a
-  C = a*izz₀
-  S₁, err, j = S₀+(e1*cⱼ-f1)*C, one(real(T)), 2
-  while err > 10eps(real(T))
-    f0, f1 = f1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*f1+z₀*(1-z₀)*(j-1)*f0+(1-2z₀)*e1+2z₀*(1-z₀)*e0)/j
-    e0, e1 = e1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*e1+z₀*(1-z₀)*(j-1)*e0)/j
-    C *= (a+j-1)*izz₀/j
-    cⱼ += 2/T(j)-one(T)/(a+j-one(T))
-    S₀, S₁ = S₁, S₁+(e1*cⱼ-f1)*C
-    err = errcheck(S₁-S₀, S₀)
-    j+=1
-  end
-  return S₁
+    T = promote_type(typeof(a), typeof(c), typeof(z₀), typeof(z))
+    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)
+    S₀ = cⱼ
+    cⱼ += 2/one(T)-one(T)/a
+    C = a*izz₀
+    S₁, j = S₀+(e1*cⱼ-f1)*C, 2
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 2
+        f0, f1 = f1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*f1+z₀*(1-z₀)*(j-1)*f0+(1-2z₀)*e1+2z₀*(1-z₀)*e0)/j
+        e0, e1 = e1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*e1+z₀*(1-z₀)*(j-1)*e0)/j
+        C *= (a+j-1)*izz₀/j
+        cⱼ += 2/T(j)-one(T)/(a+j-one(T))
+        S₀, S₁ = S₁, S₁+(e1*cⱼ-f1)*C
+        j += 1
+    end
+    return S₁
 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, S, err, j = one(T), cⱼ, one(real(T)), 0
-  while err > 10eps(real(T))
-    C *= (a+j)/(j+1)^2*(b+j)*w
-    cⱼ += 2/(j+one(T))-one(T)/(a+j)-one(T)/(b+j)
-    S += C*cⱼ
-    err = errcheck(C, S)
-    j+=1
-  end
-  return S
+    T = promote_type(typeof(a), typeof(b), typeof(z), typeof(w))
+    cⱼ = 2digamma(one(T))-digamma(a)-digamma(b)+s*log1p(-z)
+    C, S, j = one(T), cⱼ, 0
+    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
+        C *= (a+j)/(j+1)^2*(b+j)*w
+        cⱼ += 2/(j+one(T))-one(T)/(a+j)-one(T)/(b+j)
+        S += C*cⱼ
+        j += 1
+    end
+    return S
 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)
-  C, S, err, j = one(T), cⱼ, one(real(T)), 0
-  while err > 10eps(real(T))
-    C *= (a+j)/(j+1)^2*(d+j)*w
-    cⱼ += 2/(j+one(T))-one(T)/(a+j)+one(T)/(b-(j+one(T)))
-    S += C*cⱼ
-    err = errcheck(C, S)
-    j+=1
-  end
-  return S
+    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)
+    C, S, j = one(T), cⱼ, 0
+    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
+        C *= (a+j)/(j+1)^2*(d+j)*w
+        cⱼ += 2/(j+one(T))-one(T)/(a+j)+one(T)/(b-(j+one(T)))
+        S += C*cⱼ
+        j += 1
+    end
+    return S
 end
 
 function _₂F₁taylor(a::Number, b::Number, c::Number, z::Number)
-  T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
-  z₀ = abs(z) < 1 ? ρϵ*sign(z) : sign(z)/ρϵ
-  q₀, q₁ = _₂F₁(a, b, c, z₀), a*b/c*_₂F₁(a+1, b+1, c+1, z₀)
-  S₀, zz₀ = q₀, z-z₀
-  S₁, err, zz₀j, j = S₀+q₁*zz₀, one(real(T)), zz₀, 0
-  while err > 10eps(real(T))
-    q₀, q₁ = q₁, ((j*(2z₀-one(T))-c+(a+b+one(T))*z₀)*q₁ + (a+j)*(b+j)/(j+one(T))*q₀)/(z₀*(one(T)-z₀)*(j+2))
-    zz₀j *= zz₀
-    S₀, S₁ = S₁, S₁+q₁*zz₀j
-    err = errcheck(S₀-S₁, S₀)
-    j+=1
-  end
-  return S₁
+    T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
+    z₀ = abs(z) < 1 ? ρϵ*sign(z) : sign(z)/ρϵ
+    q₀, q₁ = _₂F₁(a, b, c, z₀), a*b/c*_₂F₁(a+1, b+1, c+1, z₀)
+    S₀, zz₀ = q₀, z-z₀
+    S₁, zz₀j, j = S₀+q₁*zz₀, zz₀, 0
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
+        q₀, q₁ = q₁, ((j*(2z₀-one(T))-c+(a+b+one(T))*z₀)*q₁ + (a+j)*(b+j)/(j+one(T))*q₀)/(z₀*(one(T)-z₀)*(j+2))
+        zz₀j *= zz₀
+        S₀, S₁ = S₁, S₁+q₁*zz₀j
+        j += 1
+    end
+    return S₁
 end
 
 function _₃F₂maclaurin(a₁, a₂, a₃, b₁, b₂, z)
-  T = promote_type(typeof(a₁), typeof(a₂), typeof(a₃), typeof(b₁), typeof(b₂), typeof(z))
-  S₀, S₁, err, j = one(T), one(T)+(a₁*a₂*a₃*z)/(b₁*b₂), one(real(T)), 1
-  while err > 100eps(real(T))
-    rⱼ = ((a₁+j)*(a₂+j)*(a₃+j))/((b₁+j)*(b₂+j)*(j+1))
-    S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
-    err = errcheck(S₁-S₀, S₀)
-    j+=1
-  end
-  return S₁
-end
-
-function mFnmaclaurin(a::AbstractVector{S}, b::AbstractVector{U}, z::V) where {S, U, V}
-  T = promote_type(S, U, V)
-  S₀, S₁, err, j = one(T), one(T)+prod(a)*z/prod(b), one(real(T)), 1
-  while err > 100eps(real(T))
-    rⱼ = inv(j+one(T))
-    for i=1:length(a) rⱼ *= a[i]+j end
-    for i=1:length(b) rⱼ /= b[i]+j end
-    S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
-    err = errcheck(S₁-S₀, S₀)
-    j+=1
-  end
-  return S₁
+    T = promote_type(typeof(a₁), typeof(a₂), typeof(a₃), typeof(b₁), typeof(b₂), typeof(z))
+    S₀, S₁, j = one(T), one(T)+(a₁*a₂*a₃*z)/(b₁*b₂), 1
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
+        rⱼ = ((a₁+j)*(a₂+j)*(a₃+j))/((b₁+j)*(b₂+j)*(j+1))
+        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
+        j += 1
+    end
+    return S₁
 end
 
-hypot2(a, b) = hypot(a, b)
-hypot2(a::Complex{T}, b::Complex{T}) where {T<:AbstractFloat} = hypot(hypot(a.re, b.re), hypot(a.im, b.im))
-hypot2(a::T, b::Complex{T}) where {T<:AbstractFloat} = hypot(hypot(a, b.re), b.im)
-hypot2(a::Complex{T}, b::T) where {T<:AbstractFloat} = hypot(hypot(a.re, b), a.im)
-
-errcheck(x, y) = abs(x/y)
-errcheck(x::Dual, y::Dual) = hypot2(realpart(x), dualpart(x))/hypot2(realpart(y), dualpart(y))
-
-# This is from Base because real(::Type{Dual{BigFloat}}) doesn't exist in the scop of rtoldefault, required by isapprox.
-
-# isapprox: approximate equality of numbers
-function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x, y), atol::Real=0)
-  x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y)))
+function pFqmaclaurin(a::AbstractVector{S}, b::AbstractVector{U}, z::V) where {S, U, V}
+    T = promote_type(S, U, V)
+    S₀, S₁, j = one(T), one(T)+prod(a)*z/prod(b), 1
+    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
+        rⱼ = inv(j+one(T))
+        for i=1:length(a) rⱼ *= a[i]+j end
+        for i=1:length(b) rⱼ /= b[i]+j end
+        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
+        j += 1
+    end
+    return S₁
 end
 
-const ≈ = isapprox
-≉(x, y) = !(x ≈ y)
-
-# default tolerance arguments
-rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
-rtoldefault(::Type{T}) where {T<:Real} = 0
-rtoldefault(::Type{Dual{T}}) where {T<:Real} = rtoldefault(T)
-rtoldefault(x::Union{T, Type{T}}, y::Union{S, Type{S}}) where {T<:Number, S<:Number} = rtoldefault(promote_type(real(T), real(S)))
-
-function continuedfraction(v::V, u::U, rtol::T
-    ) where {V<:Function, U<:Function, T<:Number}
-  n = 4
-  a, b = continuedfraction(v, u, n, 2n)
-  n *= 2
-  while n <= 2^16
-    isapprox(a, b, rtol=rtol) && return b
+function continuedfraction(v::V, u::U, rtol::T) where {V<:Function, U<:Function, T<:Number}
+    n = 4
+    a, b = continuedfraction(v, u, n, 2n)
     n *= 2
-    a = deepcopy(b)
-    b = continuedfraction(v, u, n)
-  end
-  isapprox(a, b, rtol=rtol) && return b
-  error("Convergence failure for continued fraction for hypergeometric function")
+    while n <= 2^16
+        isapprox(a, b, rtol=rtol) && return b
+        n *= 2
+        a = deepcopy(b)
+        b = continuedfraction(v, u, n)
+    end
+    isapprox(a, b, rtol=rtol) && return b
+    error("Convergence failure for continued fraction for hypergeometric function")
 end
-function continuedfraction(v::V, u::U, n::Int, m::Int
-    ) where {V<:Function, U<:Function}
-  @assert m > n
-  output_m = u(m) / v(m)
-  for i ∈ reverse(n+1:m-1)
-    output_m = u(i) / (v(i) + output_m)
-  end
-  ui, vi = u(n), v(n)
-  output_m = ui / (vi + output_m)
-  output_n = ui / vi
-  for i ∈ reverse(1:n-1)
-    ui, vi = u(i), v(i)
+
+function continuedfraction(v::V, u::U, n::Int, m::Int) where {V<:Function, U<:Function}
+    @assert m > n
+    output_m = u(m) / v(m)
+    for i ∈ reverse(n+1:m-1)
+        output_m = u(i) / (v(i) + output_m)
+    end
+    ui, vi = u(n), v(n)
     output_m = ui / (vi + output_m)
-    output_n = ui / (vi + output_n)
-  end
-  vi = v(0)
-  return vi + output_n, vi + output_m
+    output_n = ui / vi
+    for i ∈ reverse(1:n-1)
+        ui, vi = u(i), v(i)
+        output_m = ui / (vi + output_m)
+        output_n = ui / (vi + output_n)
+    end
+    vi = v(0)
+    return vi + output_n, vi + output_m
 end
+
 function continuedfraction(v::V, u::U, n::Int) where {V<:Function, U<:Function}
-  output = u(n) / v(n)
-  for i ∈ reverse(1:n-1)
-    output = u(i) / (v(i) + output)
-  end
-  return v(0) + output
+    output = u(n) / v(n)
+    for i ∈ reverse(1:n-1)
+        output = u(i) / (v(i) + output)
+    end
+    return v(0) + output
 end
diff --git a/src/weniger.jl b/src/weniger.jl
index 3985e0f..087ee98 100644
--- a/src/weniger.jl
+++ b/src/weniger.jl
@@ -3,9 +3,9 @@
 
 # ₘFₙ(α;β;z)
 # γ ∉ ℕ
-function wenigerpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) where {T1, T2, T3}
+function pFqweniger(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3; kmax::Int = 10_000) where {T1, T2, T3}
     T = promote_type(T1, T2, T3)
-    absα = T.(abs.(α))
+    absα = abs.(T.(α))
     if norm(z) < eps(real(T)) || norm(prod(α)) < eps(prod(absα))
         return one(T)
     end
@@ -57,7 +57,7 @@ function wenigerpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) where
     ΔD[r] = D[r+1]/pochhammer(γ-ρ-1, ρ)
     R[r+1] = N[r+1]/D[r+1]
     k = 0
-    while (abs(R[r+1]-R[r]) > 10*abs(R[r+1])*eps(real(T)) && k < 10_000) || k < r
+    while k < r || (k < kmax && errcheck(R[r], R[r+1], 10eps(real(T))))
         for j in 1:r
             N[j] = N[j+1]
             D[j] = D[j+1]
@@ -79,13 +79,13 @@ function wenigerpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) where
         end
         t2 = zero(T)
         s2 = zero(T)
-        for s in max(0,ρ+1-k):ρ
+        for s in max(0, ρ+1-k):ρ
             s2 += Cρ[s+1]*C1[s+1]*(N[r-ρ+s]+(γ+k-ρ+s-2)*N[r-ρ+s-1])
         end
         s2 += (γ+k-1)*N[r]/pochhammer(γ+2k-ρ-1, ρ+2)
         t2 += P[1]*s2
         s2 = zero(T)
-        for s in max(0,ρ+1-k):ρ+1
+        for s in max(0, ρ+1-k):ρ+1
             s2 += Cρ[s+1]*C2[s+1]*(γ+2k-2ρ+2s-3)*N[r-ρ+s-1]
         end
         ΔNold[r+1] = s2
@@ -120,12 +120,12 @@ function wenigerpFq(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3) where
         D[r+1] /= P[1]/pochhammer(γ+2k-ρ-1, ρ+2)
         R[r+1] = N[r+1]/D[r+1]
         s1 = zero(T)
-        for s in max(0,ρ-k):ρ+1
+        for s in max(0, ρ-k):ρ+1
             s1 += Cρ[s+1]*C3[s+1]*(γ+2k-2ρ+2s-1)*N[r-ρ+s]
         end
         ΔN[r+1] = s1
         s1 = zero(T)
-        for s in max(0,ρ-k):ρ+1
+        for s in max(0, ρ-k):ρ+1
             s1 += Cρ[s+1]*C3[s+1]*(γ+2k-2ρ+2s-1)*D[r-ρ+s]
         end
         ΔD[r+1] = s1
diff --git a/test/runtests.jl b/test/runtests.jl
index 08e024a..e49146b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,10 +1,10 @@
-using HypergeometricFunctions, IntervalArithmetic, Test
+using HypergeometricFunctions, Test
 import LinearAlgebra: norm
 import HypergeometricFunctions: _₂F₁general # not exported
 
 import HypergeometricFunctions: drummond0F0, drummond1F0, drummond0F1,
                                 drummond2F0, drummond1F1, drummond0F2,
-                                drummond2F1, drummondpFq, wenigerpFq
+                                drummond2F1, pFqdrummond, pFqweniger
 
 const rtol = 1.0e-3
 const NumberType = Float64
@@ -13,7 +13,7 @@ const NumberType = Float64
   @testset "_₂F₁ vs _₂F₁general" begin
     e = exp(1.0)
     regression_max_accumulated_error = 2.0e-13 # discovered by running the test
-    for z in (.9rand(Float64,10), 10rand(ComplexF64, 10))
+    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),
                       (-0.3, 1.3, 0.5), (0.35, 0.85, 0.5), (0.5, 0.5, 1.5),
@@ -36,188 +36,188 @@ const NumberType = Float64
   end
 
 
-  @testset "mFn vs mpmath" begin
+  @testset "pFq vs mpmath" begin
     (a, b, c, result) = NumberType.([-1.138]), NumberType.([-1.17865]), NumberType(0.29524203533943405), 1.3211939223293958
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.00209496]), NumberType.([1.55444]), NumberType(1.9964304489951332), 0.9956895798837077
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.550861]), NumberType.([-2.76028]), NumberType(1.686164802648714), 2.26529646755787
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.522481]), NumberType.([2.43767]), NumberType(0.07397564666157663), 1.0161190845825492
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.14629]), NumberType.([-1.71237]), NumberType(0.865999806211927), 6.09627259912178
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.977563]), NumberType.([-0.943489]), NumberType(1.3295519796473885), 2.9650972591473237
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.19487]), NumberType.([-2.84673]), NumberType(-2.946765152780697), -10.159918912696533
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.40648]), NumberType.([-2.09642]), NumberType(0.003973775403141477), 1.002667944757093
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.78115]), NumberType.([-0.974578]), NumberType(-0.289876068484225), -9.436220915244991
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.0195]), NumberType.([2.97415]), NumberType(1.9088632413219395), 0.020342030150391748
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.29058]), NumberType.([1.12475, 2.20863]), NumberType(2.5233777332943026), -0.23709345988925976
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.93422]), NumberType.([-0.403018, -1.68818]), NumberType(1.8841553851312423), -368.4667336251361
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.38965]), NumberType.([1.34519, 2.18274]), NumberType(-2.967668498991821), 0.12382284309201524
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.317445]), NumberType.([2.83319, 2.08577]), NumberType(-0.8292545561139764), 0.9574454648148745
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.53034]), NumberType.([-2.00218, 2.71893]), NumberType(-0.055950986916451395), 0.984386941287138
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.8068]), NumberType.([2.785, -2.27362]), NumberType(-2.965224079697172), 3.343788228467509
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.01277]), NumberType.([0.0297853, 0.638178]), NumberType(2.87029742430684), 1862.664910896795
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.332]), NumberType.([-0.571085, 2.40734]), NumberType(2.5191976052414615), 2.6748631921200974
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.585852]), NumberType.([-0.840541, -0.914398]), NumberType(-0.7114454462840016), 13.397306612348904
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.30148]), NumberType.([-0.0397773, 0.452281]), NumberType(2.9085746526416285), 134.70529132344416
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.79664]), NumberType.([2.5454, -2.0446, -1.84776]), NumberType(-1.5647341841975129), 10.224086634371103
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.08954]), NumberType.([-1.09992, -0.0898055, -1.21452]), NumberType(2.911146030454528), -281.8950930467646
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.50123]), NumberType.([1.481, 1.04535, -1.1706]), NumberType(1.2093267908720815), 4.297922628682247
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.918156]), NumberType.([-1.66879, -1.26548, -2.71329]), NumberType(-1.4539176485483885), 0.7245506668664914
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.25871]), NumberType.([0.734527, 2.69253, 1.52627]), NumberType(1.412857816927049), 0.0007367729643306418
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.75347]), NumberType.([-1.19881, 1.46644, 0.377671]), NumberType(-2.894612064681616), 31.986877557865302
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.55356]), NumberType.([2.31273, -1.79583, 2.71976]), NumberType(-1.8211352681511421), 0.6524733226915893
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.96471]), NumberType.([0.10731, 0.150247, -0.749495]), NumberType(-0.5456496601437784), -364.7310910695041
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.28595]), NumberType.([1.24272, -1.3733, -1.97931]), NumberType(2.694315046446391), 1427.85053117312
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.02891]), NumberType.([-0.203779, 0.414445, -2.72581]), NumberType(2.7264202122777323), -11.325191332268156
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.21643, 2.56177]), NumberType.([-2.85355]), NumberType(-1.173143685351698), -2.6927101990763944
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.884016, -0.278806]), NumberType.([-2.95573]), NumberType(-0.5131107181888153), 1.0423537583061855
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.489218, 1.02043]), NumberType.([-1.8805]), NumberType(-2.054787166536263), 1.3985297910285344
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.61829, -0.0405846]), NumberType.([-2.52321]), NumberType(-2.8547909742932127), 1.111640235796153
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.78806, 0.244808]), NumberType.([-2.03175]), NumberType(0.5005127629084343), 1.0071139633716848
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.79663, 2.95841]), NumberType.([-2.23659]), NumberType(-2.1658329294772676), -1.298805755547899
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.473781, -0.844289]), NumberType.([-0.228894]), NumberType(0.8762669945656332), 3.026596343521232
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.99653, -2.15667]), NumberType.([1.45492]), NumberType(0.5669144675224369), 2.901074111158336
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.89145, -2.17305]), NumberType.([-1.97603]), NumberType(-2.2566276457677343), 145.35756041830106
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.998819, -0.225309]), NumberType.([0.447173]), NumberType(-0.10835817353459554), 1.051599100560819
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.383828, 1.77538]), NumberType.([-1.0875, -2.84315]), NumberType(1.4065086054194986), 715.336897395762
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.733652, -1.6307]), NumberType.([-0.0608787, 0.896927]), NumberType(0.10210464529041374), 3.1665794151940694
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.202523, -1.68826]), NumberType.([-1.92541, 2.88428]), NumberType(0.23168614635513318), 0.9854522687722247
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.946942, 1.58094]), NumberType.([-0.150083, 2.30757]), NumberType(-2.109785441324441), -7.76796376651872
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.55118, -2.75387]), NumberType.([1.46214, 2.18634]), NumberType(2.6348816640393258), -0.31994633361734
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.06745, -0.28521]), NumberType.([-1.99816, -0.766774]), NumberType(-0.7965581191453945), 5.472490531872681
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.20859, -1.41744]), NumberType.([1.85829, 1.92196]), NumberType(-1.2487286295352975), -0.053122867013287486
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.260072, 0.317936]), NumberType.([-1.4299, 1.60909]), NumberType(-1.554201646172852), 1.0698093358492726
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.81949, 2.65844]), NumberType.([2.84994, -1.67876]), NumberType(1.5573520464968), 4.557008143231845
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.09665, 1.27147]), NumberType.([-0.682723, 0.192544]), NumberType(-1.175497640118377), -118.5195721586634
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.472347, 2.38973]), NumberType.([1.23953, -2.45145, -0.0816949]), NumberType(2.9793981957869673), 562.1327877798817
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.70886, -2.26978]), NumberType.([-1.83008, -1.6277, -1.33282]), NumberType(2.7790611310114053), -12231.7908950255
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.48786, -2.73106]), NumberType.([-1.73601, -1.676, -2.17108]), NumberType(-0.6847891664895385), -36.681955083411964
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.72604, -1.83087]), NumberType.([-2.68909, -0.904331, 1.93759]), NumberType(1.651593934432523), 0.9572670911067562
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.27953, 2.57139]), NumberType.([-1.72806, 2.13894, -0.674948]), NumberType(1.3654741919794335), -8.061422671050726
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.363687, 1.0458]), NumberType.([-1.13771, -2.02746, -2.86789]), NumberType(2.026179003725624), -5919.391531381284
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.49209, -1.06563]), NumberType.([1.47807, 1.01142, -1.54966]), NumberType(-0.9883650134564901), -0.0985713724396932
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.416307, 1.24566]), NumberType.([1.56196, 1.32576, -1.82953]), NumberType(-2.4172898932733053), 1.0404646069316064
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.52485, 0.109503]), NumberType.([-1.3165, 0.177332, 2.58557]), NumberType(-2.84982197723157), 2.9369917768818667
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.45525, -1.61598]), NumberType.([-1.75953, -1.59772, 2.05578]), NumberType(-1.1643127384783538), 0.5520841179227342
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.7029, 1.4496, -2.82896]), NumberType.([-2.77533]), NumberType(-1.5171018909198049), 8.437330705993508
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.34783, -0.324078, 2.44238]), NumberType.([-0.777003]), NumberType(-1.059855902026896), 3.9754790978546266
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.87289, 2.99671, 2.99742]), NumberType.([-1.4312]), NumberType(-2.433424057786612), 638.863801591021
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.963631, -1.36196, -2.90998]), NumberType.([-0.950805]), NumberType(-1.1275151312466227), -1.9506846998048282
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.20585, 1.21729, 0.419755]), NumberType.([-2.95176]), NumberType(-0.8474785819006732), -0.9897576725972367
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.50248, -2.01281, -1.10468]), NumberType.([-1.47035]), NumberType(-2.161852515937607), -4.205644305802463
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.373, 0.0281506, 1.31409]), NumberType.([-0.775952]), NumberType(-2.0586457735217407), 0.6176889810366223
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.4476, 0.286052, 1.71704]), NumberType.([1.29826]), NumberType(-0.9902327227931127), 0.7640485766337997
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.49284, 1.34398, 1.38886]), NumberType.([-2.04354]), NumberType(-0.9516109943455211), -137.34441283528875
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.19033, -0.446436, -2.73578]), NumberType.([0.758367]), NumberType(-1.115415216535085), -8.890490669826645
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.21084, -1.43893, 1.56842]), NumberType.([-0.598469, 0.52221]), NumberType(-2.2391561626302376), -64.66480492684083
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.607277, -0.376803, -2.08424]), NumberType.([1.92106, 2.80899]), NumberType(-1.5589765316931263), 0.8516398920401931
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.63828, 2.59771, 2.54461]), NumberType.([2.62021, -0.0157026]), NumberType(0.5526474385388149), -29.571181883739737
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.314356, 2.59205, -1.28807]), NumberType.([-1.22202, -1.66475]), NumberType(-0.7121680260550485), 0.856714022601912
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.694337, -2.06707, 0.214575]), NumberType.([-1.31576, 0.266519]), NumberType(-2.9461185849069706), 21.80248076992691
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.90394, 2.56441, 0.07523]), NumberType.([0.0390543, -1.65175]), NumberType(-0.7780439739130411), -3.184516491992518
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-2.80752, 2.7876, -1.37287]), NumberType.([-2.9165, 0.49188]), NumberType(0.39686964583141693), -1.3014023309227971
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([2.36231, -0.901808, 1.00893]), NumberType.([1.39985, -2.34968]), NumberType(0.970806084927347), 42731.87930657499
-    @test_broken mFn(a, b, c) ≈ result atol=eps() rtol=rtol
-    @test HypergeometricFunctions.mFncontinuedfraction(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test_broken pFq(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test HypergeometricFunctions.pFqcontinuedfraction(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.715302, 2.94101, -0.498188]), NumberType.([-0.393808, -0.0265878]), NumberType(0.4375358657735875), 61.68955383799459
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.14105, 1.9464, -2.44787]), NumberType.([0.54636, -0.654353]), NumberType(-2.9505198429994763), -53.189709409599516
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([0.228056, -1.24092, 0.843438]), NumberType.([1.25348, 2.42447, -1.59626]), NumberType(-2.32373843841519), 0.8942755488371047
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.78019, -2.49326, -1.62843]), NumberType.([-1.85186, -2.72869, -2.36236]), NumberType(-1.6400906330598044), 0.2981496993840609
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.98727, 2.31523, -0.0735596]), NumberType.([-0.0644426, 2.54497, -1.08036]), NumberType(-1.0705461999772998), 10.575143699481718
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.45036, -2.25443, -2.74298]), NumberType.([-1.39353, -1.20841, 0.225826]), NumberType(-2.950872722868197), 982.2060223576477
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-0.36497, 1.79508, 2.36814]), NumberType.([1.55155, 0.704328, -0.686943]), NumberType(0.2186808938362801), 1.7064306875530209
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.51876, -1.37851, -1.12183]), NumberType.([1.3349, -2.95129, 1.10139]), NumberType(0.7674847040296018), 1.4158192973160308
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.99362, -2.09819, -2.16471]), NumberType.([-1.76093, 2.0683, -0.947339]), NumberType(-1.8729057277220629), -41.660320326270636
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.82466, 1.28428, 2.67447]), NumberType.([2.10754, 0.871646, 1.39154]), NumberType(2.658591557834844), 81.52661900335721
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([-1.79631, -0.0269564, 0.915417]), NumberType.([-2.17853, 1.0933, -0.000243005]), NumberType(2.661306310651646), -83.69268822607185
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
     (a, b, c, result) = NumberType.([1.36163, -0.758606, -0.716056]), NumberType.([-1.04856, 0.266574, -2.47814]), NumberType(1.8461964795702661), -113.75411958323524
-    @test mFn(a, b, c) ≈ result atol=eps() rtol=rtol
+    @test pFq(a, b, c) ≈ result atol=eps() rtol=rtol
   end
   @testset "_₂F₁ vs mpmath" begin
     (a, b, c, result) = NumberType.([-0.509718, 2.74761]), NumberType.([-0.650532]), NumberType(-2.9555898601070316), 1.2790428246176058
@@ -271,11 +271,21 @@ end
                    (Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for z in S(-4):S(0.25):S(1)
-            @test drummond0F0(z) ≈ S(mFn(T[], T[], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[], S[], z) ≈ S(mFn(T[], T[], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[], S[], z) ≈ S(mFn(T[], T[], T(z))) atol=atol rtol=rtol
+            @test drummond0F0(z) ≈ S(pFq(T[], T[], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[], z) ≈ S(pFq(T[], T[], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[], z) ≈ S(pFq(T[], T[], T(z))) atol=atol rtol=rtol
+        end
+        atol *= 2
+        rtol *= 2
+        for x in S(-4):S(0.25):S(1), y in S(-4):S(0.25):S(1)
+            z = complex(x, y)
+            @test drummond0F0(z) ≈ CS(pFq(T[], T[], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[], z) ≈ CS(pFq(T[], T[], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[], z) ≈ CS(pFq(T[], T[], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -284,11 +294,21 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for α in S(-1.5):S(0.5):S(1.5), z in S(-0.75):S(0.25):S(0.25)
-            @test drummond1F0(α, z) ≈ S(mFn(T[α], T[], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[α], S[], z) ≈ S(mFn(T[α], T[], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[α], S[], z) ≈ S(mFn(T[α], T[], T(z))) atol=atol rtol=rtol
+            @test drummond1F0(α, z) ≈ S(pFq(T[α], T[], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α], S[], z) ≈ S(pFq(T[α], T[], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α], S[], z) ≈ S(pFq(T[α], T[], T(z))) atol=atol rtol=rtol
+        end
+        atol *= 2
+        rtol *= 2
+        for α in S(-1.5):S(0.5):S(1.5), x in S(-0.75):S(0.25):S(0.25), y in S(-0.75):S(0.25):S(0.25)
+            z = complex(x, y)
+            @test drummond1F0(α, z) ≈ CS(pFq(T[α], T[], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α], S[], z) ≈ CS(pFq(T[α], T[], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α], S[], z) ≈ CS(pFq(T[α], T[], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -297,11 +317,21 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for α in S(-1.5):S(1.0):S(1.5), z in S(-0.75):S(0.25):S(0.75)
-            @test drummond0F1(α, z) ≈ S(mFn(T[], T[α], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[], S[α], z) ≈ S(mFn(T[], T[α], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[], S[α], z) ≈ S(mFn(T[], T[α], T(z))) atol=atol rtol=rtol
+            @test drummond0F1(α, z) ≈ S(pFq(T[], T[α], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[α], z) ≈ S(pFq(T[], T[α], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[α], z) ≈ S(pFq(T[], T[α], T(z))) atol=atol rtol=rtol
+        end
+        atol *= 2
+        rtol *= 2
+        for α in S(-1.5):S(1.0):S(1.5), x in S(-0.75):S(0.25):S(0.75), y in S(-0.75):S(0.25):S(0.75)
+            z = complex(x, y)
+            @test drummond0F1(α, z) ≈ CS(pFq(T[], T[α], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[α], z) ≈ CS(pFq(T[], T[α], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[α], z) ≈ CS(pFq(T[], T[α], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -310,11 +340,19 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
-        atol = rtol = 1000eps(S)
+        CS = Complex{S}
+        CT = Complex{T}
+        atol = rtol = sqrt(eps(S))
         for α in S(-1.5):S(1.0):S(1.5), β in S(-1.5):S(1.0):S(1.5), z in S(-1.0):S(0.25):S(0.0)
             @test drummond2F0(α, β, z) ≈ S(drummond2F0(T(α), T(β), T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[α, β], S[], z) ≈ S(drummondpFq(T[α, β], T[], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[α, β], S[], z) ≈ S(wenigerpFq(T[α, β], T[], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α, β], S[], z) ≈ S(pFqdrummond(T[α, β], T[], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α, β], S[], z) ≈ S(pFqweniger(T[α, β], T[], T(z))) atol=atol rtol=rtol
+        end
+        for α in S(-0.5):S(1.0):S(0.5), β in S(-0.5):S(1.0):S(0.5), x in S(-0.5):S(0.25):S(0.0), y in S(-0.5):S(0.25):S(0.0)
+            z = complex(x, y)
+            @test drummond2F0(α, β, z) ≈ CS(drummond2F0(T(α), T(β), CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α, β], S[], z) ≈ CS(pFqdrummond(T[α, β], T[], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α, β], S[], z) ≈ CS(pFqweniger(T[α, β], T[], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -323,11 +361,21 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for α in S(-1.5):S(0.5):S(1.5), β in S(-1.5):S(1.0):S(1.5), z in S(-0.75):S(0.25):S(0.75)
-            @test drummond1F1(α, β, z) ≈ S(mFn(T[α], T[β], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[α], S[β], z) ≈ S(mFn(T[α], T[β], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[α], S[β], z) ≈ S(mFn(T[α], T[β], T(z))) atol=atol rtol=rtol
+            @test drummond1F1(α, β, z) ≈ S(pFq(T[α], T[β], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α], S[β], z) ≈ S(pFq(T[α], T[β], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α], S[β], z) ≈ S(pFq(T[α], T[β], T(z))) atol=atol rtol=rtol
+        end
+        atol *= 2
+        rtol *= 2
+        for α in S(-1.5):S(0.5):S(1.5), β in S(-1.5):S(1.0):S(1.5), x in S(-0.75):S(0.25):S(0.75), y in S(-0.75):S(0.25):S(0.75)
+            z = complex(x, y)
+            @test drummond1F1(α, β, z) ≈ CS(pFq(T[α], T[β], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α], S[β], z) ≈ CS(pFq(T[α], T[β], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α], S[β], z) ≈ CS(pFq(T[α], T[β], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -336,11 +384,21 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for α in S(-1.5):S(1.0):S(1.5), β in S(-1.5):S(1.0):S(1.5), z in S(-0.75):S(0.25):S(0.75)
-            @test drummond0F2(α, β, z) ≈ S(mFn(T[], T[α, β], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[], S[α, β], z) ≈ S(mFn(T[], T[α, β], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[], S[α, β], z) ≈ S(mFn(T[], T[α, β], T(z))) atol=atol rtol=rtol
+            @test drummond0F2(α, β, z) ≈ S(pFq(T[], T[α, β], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[α, β], z) ≈ S(pFq(T[], T[α, β], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[α, β], z) ≈ S(pFq(T[], T[α, β], T(z))) atol=atol rtol=rtol
+        end
+        atol *= 2
+        rtol *= 2
+        for α in S(-1.5):S(1.0):S(1.5), β in S(-1.5):S(1.0):S(1.5), x in S(-0.75):S(0.25):S(0.75), y in S(-0.75):S(0.25):S(0.75)
+            z = complex(x, y)
+            @test drummond0F2(α, β, z) ≈ CS(pFq(T[], T[α, β], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[], S[α, β], z) ≈ CS(pFq(T[], T[α, β], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[], S[α, β], z) ≈ CS(pFq(T[], T[α, β], CT(z))) atol=atol rtol=rtol
         end
     end
 end
@@ -349,26 +407,21 @@ end
     for (S, T) in ((Float32, Float64),
                    (Float64, BigFloat),
                    (BigFloat, BigFloat))
+        CS = Complex{S}
+        CT = Complex{T}
         atol = rtol = 1000eps(S)
         for α in S(-1.5):S(0.5):S(1.5), β in S(-1.5):S(0.5):S(1.5), γ in S(-1.5):S(1.0):S(1.5), z in S(-0.625):S(0.25):S(0.125)
-            @test drummond2F1(α, β, γ, z) ≈ S(mFn(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
-            @test drummondpFq(S[α, β], S[γ], z) ≈ S(mFn(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
-            @test wenigerpFq(S[α, β], S[γ], z) ≈ S(mFn(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
+            @test drummond2F1(α, β, γ, z) ≈ S(pFq(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α, β], S[γ], z) ≈ S(pFq(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α, β], S[γ], z) ≈ S(pFq(T[α, β], T[γ], T(z))) atol=atol rtol=rtol
         end
-    end
-end
-
-@testset "₂F₁ with interval arithmetic" begin
-    for (S, IS, T, IT) in ((Float32, Interval{Float32}, Float64, Interval{Float64}),
-                   (Float64, Interval{Float64}, BigFloat, Interval{BigFloat}),
-                   (BigFloat, Interval{BigFloat}, BigFloat, Interval{BigFloat}))
-        atol = rtol = sqrt(eps(S))
-        α = S(0.25)
-        β = S(0.5)
-        γ = S(0.375)
-        for z in S(-0.125):S(0.25):S(0.125)
-            @test drummondpFq(IS[α, β], IS[γ], IS(z)) ≈ IS(mFn(IT[α, β], IT[γ], IT(z))) atol=atol rtol=rtol
-            @test wenigerpFq(IS[α, β], IS[γ], IS(z)) ≈ IS(mFn(IT[α, β], IT[γ], IT(z))) atol=atol rtol=rtol
+        atol *= 2
+        rtol *= 2
+        for α in S(-1.5):S(0.5):S(1.5), β in S(-1.5):S(0.5):S(1.5), γ in S(-1.5):S(1.0):S(1.5), x in S(-0.375):S(0.25):S(0.125), y in S(-0.375):S(0.25):S(0.125)
+            z = complex(x, y)
+            @test drummond2F1(α, β, γ, z) ≈ CS(pFq(T[α, β], T[γ], CT(z))) atol=atol rtol=rtol
+            @test pFqdrummond(S[α, β], S[γ], z) ≈ CS(pFq(T[α, β], T[γ], CT(z))) atol=atol rtol=rtol
+            @test pFqweniger(S[α, β], S[γ], z) ≈ CS(pFq(T[α, β], T[γ], CT(z))) atol=atol rtol=rtol
         end
     end
 end