fix common parameter cancelation

Project.toml

@@ -1,6 +1,6 @@
 name = "HypergeometricFunctions"
 uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
-version = "0.3.29"
+version = "0.3.30"
 
 [deps]
 Gamma = "a0844989-3bd2-4988-8bea-c9407ab0941b"

src/conformal.jl

@@ -283,7 +283,7 @@ function pFqconformalrational(α::Tuple{T1, Vararg{T1, q}}, β::NTuple{q, T2}, z
     α = Tuple(2*m .* T.(α))
     β = Tuple(2*m .* (T.(β) .- 1))
 
-    # Defining factors L_k and R_k as operators     # L_0...L_q w= R_1...R_{q+1}
+    # Defining factors L_k and R_k as operators     # L_0⋯L_q w = R_1⋯R_{q+1} w
     L = zeros(T, m+2, q)
     for j in 1:q-1
         for k in 0:m+1
@@ -415,3 +415,141 @@ function pFqconformalrational(α::Tuple{T1, Vararg{T1, q}}, β::NTuple{q, T2}, z
 
     return Shi
 end
+
+"""
+Maclaurin series of ₚ₊₁Fₚ(α, β; z) via the generalized Euler transformation with γ = ∑β-∑α:
+
+ₚ₊₁Fₚ(α, β; z) = (1-z)^γ w(z = 1 - ((1-ζ)/(1+ζ))^m)
+
+"""
+pFqconformalrationalgeneralizedeuler(α, β, z; kwds...) = pFqconformalrationalgeneralizedeuler(α, β, z, Val(4); kwds...)
+function pFqconformalrationalgeneralizedeuler(α::Tuple{Any, Vararg{Any, q}}, β::NTuple{q, Any}, z, m; kwds...) where q
+    T1 = isempty(α) ? Any : mapreduce(typeof, promote_type, α)
+    T2 = isempty(β) ? Any : mapreduce(typeof, promote_type, β)
+    pFqconformalrationalgeneralizedeuler(T1.(α), T2.(β), z, m; kwds...)
+end
+
+
+function pFqconformalrationalgeneralizedeuler(α::Tuple{T1, Vararg{T1, q}}, β::NTuple{q, T2}, z::T3, ::Val{m}; kmax::Int = KMAX) where {q, m, T1, T2, T3}
+    T = promote_type(T1, T2, T3)
+    if abs(z) < eps(real(T))
+        return one(T)
+    end
+    γ = sum(β)-sum(α)
+    # Transform alpha and beta
+    β = Tuple(2*m .* (T.(β) .- 1))
+    α = Tuple(2*m .* T.(α))
+    # Defining factors L_k and R_k as operators     # L_0⋯L_q w = R_1⋯R_{q+1} w
+
+    L = zeros(T, m+2, q)
+    for j in 1:q
+        for k in 0:m+1
+            L[k+1, j] = (binomial(m, k+1) * (1 - (-1)^(k+1)) - binomial(m, k-1) * (1 - (-1)^(k-1))) * (0 - k)- (1 - (-1)^k)* binomial(m, k) *((q - j) * m + β[j]+2m*(γ-q+j))+ (1 - (-1)^(k-1))*binomial(m, k-1)*(q-j)*m+ β[j] * binomial(m, k)
+        end
+    end
+    ΔL = zeros(T, m+2)
+    for k in 0:2:m+1
+        ΔL[k+1] = 2*(binomial(m, k+1) - binomial(m, k-1))
+    end
+
+    R = zeros(T, m+2, q+1)
+    for j in 1:q+1
+        for k in 0:m+1
+            R[k+1, j] = (binomial(m, k+1) * (1 - (-1)^(k+1)) - binomial(m, k-1) * (1 - (-1)^(k-1))) * (0 - k) - (1 - (-1)^k)*binomial(m, k) *((q - j + 1) *m + α[j] + 2m*(γ-q+j-1)) +  (1 - (-1)^(k-1))*binomial(m, k-1)*(q-j+1)*m+ α[j] * binomial(m, k)
+        end
+    end
+    ΔR = zeros(T, m+2)
+    for k in 0:2:m+1
+        ΔR[k+1] = 2*(binomial(m, k+1) - binomial(m, k-1))
+    end
+
+    # Truncated L_0
+    Lt0 = zeros(T, m+3)
+    for k in 0:m+2
+        Lt0[k+1] = (binomial(m, k+1) - binomial(m, k-1)) * (0 - k) + q*m*binomial(m, k-1)+m*(q-2γ)*binomial(m,k)
+    end
+    ΔLt0 = zeros(T, m+2)
+    for k in 0:m+1
+        ΔLt0[k+1] = binomial(m, k+1) - binomial(m, k-1)
+    end
+
+    # Initial values for the sum
+    ζ   = r(T(z), Val(m))
+    Slo = zero(T)
+    Shi = one(T)
+
+    ##############################################################################################################
+    # Matrices of auxiliary sequences, Z and Zt
+    Z  = zeros(T, m+2, q+1)
+    Zt = zeros(T, m+2, q)
+    # Initialize Z (n=0)
+    Z[m+2, 1] = one(T)
+    for i in 0:q-1
+        Z[m+2, i+2] = sum(L[k+1, q-i]*Z[m+3-(k+1), i+1] for k in 0:m+1)/ζ
+    end
+
+    # Initialize Zt (n=0)
+    Zt[m+2, 1] = sum(R[k+1, q+1]*Z[m+3-(k+1), 1] for k in 0:m+1)
+    for i in 1:q-1
+        Zt[m+2, i+1] = sum(R[k+1, q+1-i]*Zt[m+3-(k+1), i] for k in 0:m+1)/ζ
+    end
+
+    #############   Main loop updating matrices and summing
+
+    n = 0
+    @inbounds while (n < kmax && errcheck(Shi, Slo, 8eps(real(T))))
+        # Increment L
+        for j in 1:q
+            for k in 0:2:m+1
+                L[k+1, j] += ΔL[k+1]
+            end
+        end
+
+        # Increment Z
+        vn_new = (sum(R[k+1, 1]*Zt[m+3-(k+1), q] for k in 0:m+1) - ζ*sum(Lt0[k+1]*Z[m+3-(k+1), q+1] for k in 0:m+1)) / (n+1)
+        for j in 1:m+1
+            Z[j, q+1] = ζ*Z[j+1, q+1]
+        end
+        Z[m+2, q+1] = vn_new
+
+        for i in 1:q
+            zni_new = ζ*(Z[m+2, q+1-i+1] - sum(L[k+1, i]*Z[m+4-(k+1), q+1-i] for k in 1:m+1)) / (2*(n+1)*m + β[i])
+            for j in 1:m+1
+                Z[j, q+1-i] = ζ*Z[j+1, q+1-i]
+            end
+            Z[m+2, q+1-i] = zni_new
+        end
+
+        # Increment Lt0
+        for k in 0:m+1
+            Lt0[k+1] += ΔLt0[k+1]
+        end
+
+        # Increment R
+        for j in 1:q+1
+            for k in 0:2:m+1
+                R[k+1, j] += ΔR[k+1]
+            end
+        end
+
+        # Increment Zt
+        zt1_new = sum(R[k+1, q+1]*Z[m+3-(k+1), 1] for k in 0:m+1)
+        for j in 1:m+1
+            Zt[j, 1] = ζ*Zt[j+1, 1]
+        end
+        Zt[m+2, 1] = zt1_new
+
+        for i in 1:q-1
+            ztni_new = sum(R[k+1, q+1-i]*Zt[m+3-(k+1), i] for k in 0:m+1)/ζ
+            for j in 1:m+1
+                Zt[j, i+1] = ζ*Zt[j+1, i+1]
+            end
+            Zt[m+2, i+1] = ztni_new
+        end
+
+        Shi, Slo = Shi + Z[m+2, 1], Shi
+        n += 1
+    end
+
+    return Shi*exp(γ*log1p(-z))
+end

src/gauss.jl

@@ -84,7 +84,7 @@ function _₂F₁general(a, b, c, z; kwds...)
         return _₂F₁general(b, a, c, z; kwds...)
     elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
         return _₂F₁maclaurin(a, b, c, z; kwds...)
-    elseif !isalmostwellpoised(a, b, c)
+    elseif ispositivelybalanced(a, b, c)
         return exp((c-a-b)*log1p(-z))*_₂F₁general(c-a, c-b, c, z; kwds...)
     elseif abs(z/(z-1)) ≤ ρ
         return exp(-a*log1p(-z))_₂F₁maclaurin(a, c-b, c, z/(z-1); kwds...)
@@ -115,7 +115,7 @@ function _₂F₁general2(a, b, c, z; kwds...)
         return _₂F₁general2(b, a, c, z; kwds...)
     elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀
         return _₂F₁maclaurin(a, b, c, z; kwds...)
-    elseif !isalmostwellpoised(a, b, c)
+    elseif ispositivelybalanced(a, b, c)
         return exp((c-a-b)*log1p(-z))*_₂F₁general2(c-a, c-b, c, z; kwds...)
     elseif abs(z / (z - 1)) ≤ ρ && absarg(1 - z) < convert(real(T), π) # 15.8.1
         w = z/(z-1)
@@ -175,9 +175,9 @@ end
 function _₂F₁argument_unity(a, b, c, z; kwds...)
     T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z))
     a, b, c = T(a), T(b), T(c)
-    if iswellpoised(a, b, c)
+    if ispositivelybalanced(a, b, c)
         return exp(loggamma(c)-loggamma(c-a)+loggamma(c-a-b)-loggamma(c-b))
-    elseif isalmostwellpoised(a, b, c)
+    elseif isalmostpositivelybalanced(a, b, c)
         if a + b == c
             f = loggamma(c)-loggamma(a)-loggamma(b)
             if isreal(f)
@@ -188,7 +188,7 @@ function _₂F₁argument_unity(a, b, c, z; kwds...)
         else
             return T(NaN)
         end
-    else # !isalmostwellpoised(a, b, c)
+    else # !isalmostpositivelybalanced(a, b, c)
         f = loggamma(c)-loggamma(a)+gamma(a+b-c)-loggamma(b)
         if isreal(f)
             return T(Inf)

src/generalized.jl

@@ -58,69 +58,58 @@ pFq(α::NTuple{1, <: Integer}, ::Tuple{}, z; kwds...) = (1-z)^-α[1]
 pFq(α::NTuple{1}, β::NTuple{1}, z; kwds...) = _₁F₁(α[1], β[1], z; kwds...)
 pFq(α::NTuple{2, Any}, β::NTuple{1}, z; kwds...) = _₂F₁(α[1], α[2], β[1], z; kwds...)
 
-"""
-Remove parameters common to `α` and `β`, since `(x)ₖ/(x)ₖ ≡ 1` for any shared `x`.
-Cancelling them avoids loss of accuracy (or a `0/0`,
-if `x` is a nonpositive integer) in the generic algorithms.
-"""
-cancelcommonparameters(α::Tuple{}, β::Tuple) = (α, β)
-function cancelcommonparameters(α::Tuple, β::Tuple)
-    a, resta = first(α), Base.tail(α)
-    β′ = _deletefirstequal(a, β)
-    if length(β′) == length(β)
-        resta′, β″ = cancelcommonparameters(resta, β)
-        return (a, resta′...), β″
+function pFq(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q}
+    T1 = isempty(α) ? Any : mapreduce(typeof, promote_type, α)
+    T2 = isempty(β) ? Any : mapreduce(typeof, promote_type, β)
+    pFq(T1.(α), T2.(β), z; kwds...)
+end
+function pFq(α::NTuple{p, T1}, β::Tuple{}, z; kwds...) where {p, T1}
+    z = float(z)
+    T = promote_type(T1, typeof(z))
+    if isterminating(α)
+        return pFqmaclaurin(α, β, z, getterminatinglength(α); kwds...)::T
     else
-        return cancelcommonparameters(resta, β′)
+        return pFqweniger(α, β, z; kwds...)::T
     end
 end
-
-"""
-Remove the first element of `t` that equals `x`, or return `t` unchanged if none does.
-"""
-_deletefirstequal(x, t::Tuple{}) = t
-function _deletefirstequal(x, t::Tuple)
-    y, rest = first(t), Base.tail(t)
-    return y == x ? rest : (y, _deletefirstequal(x, rest)...)
-end
-
-function pFq(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q}
-    α′, β′ = cancelcommonparameters(α, β)
-    if length(α′) != p || length(β′) != q
-        return pFq(α′, β′, z; kwds...)
+function pFq(α::Tuple{}, β::NTuple{q, T2}, z; kwds...) where {q, T2}
+    z = float(z)
+    T = promote_type(T2, typeof(z))
+    if real(z) > 0
+        return pFqmaclaurin(α, β, z; kwds...)::T
+    else
+        return pFqweniger(α, β, z; kwds...)::T
     end
+end
+function pFq(α::NTuple{p, T1}, β::NTuple{q, T2}, z; kwds...) where {p, q, T1, T2}
     z = float(z)
+    T = promote_type(T1, T2, typeof(z))
+    if isacommonparameter(α, β)
+        return pFq(cancelacommonparameter(α, β)..., z; kwds...)::T
+    end
+    if isterminating(α)
+        return pFqmaclaurin(α, β, z, getterminatinglength(α); kwds...)::T
+    end
     if p ≤ q
         if real(z) > 0
-            return pFqmaclaurin(α, β, z; kwds...)
+            return pFqmaclaurin(α, β, z; kwds...)::T
         else
-            return pFqweniger(α, β, z; kwds...)
+            return pFqweniger(α, β, z; kwds...)::T
         end
     elseif p == q + 1
-        if abs(z) ≤ ρ
-            return pFqmaclaurin(α, β, z; kwds...)
+        if ispositivelybalanced(α, β)
+            #return pFqconformalrationalgeneralizedeuler(α, β, z; kwds...)::T
+            return pFqconformalrational(α, β, z; kwds...)::T
         else
-            return pFqweniger(α, β, z; kwds...)
+            return pFqconformalrational(α, β, z; kwds...)::T
         end
     else
-        return pFqweniger(α, β, z; kwds...)
+        return pFqweniger(α, β, z; kwds...)::T
     end
 end
 
 pFq(α::AbstractVector, β::AbstractVector, z; kwds...) = pFq(Tuple(α), Tuple(β), z; kwds...)
 
-"""
-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₁, 1, 1, b₁, 2, z)`.
 """