From e4d40eda92a2a025f4c197d97fcda6d56a2f96ea Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Wed, 29 Jul 2026 11:43:29 -0500 Subject: [PATCH 1/3] fix common parameter cancelation add terminating check add generalized euler transformation (in conformal rational approximation) --- src/conformal.jl | 140 +++++++++++++++++++++++++++++++++++++++- src/gauss.jl | 10 +-- src/generalized.jl | 127 ++++++++++++++++++++++++++++++++---- src/specialfunctions.jl | 96 ++++++++++++++++++++++++--- test/runtests.jl | 64 +++++++++++------- 5 files changed, 387 insertions(+), 50 deletions(-) diff --git a/src/conformal.jl b/src/conformal.jl index 04e3db5..8bd9b07 100644 --- a/src/conformal.jl +++ b/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 diff --git a/src/gauss.jl b/src/gauss.jl index a0fea29..817d717 100644 --- a/src/gauss.jl +++ b/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) diff --git a/src/generalized.jl b/src/generalized.jl index 12b975e..70aba3b 100644 --- a/src/generalized.jl +++ b/src/generalized.jl @@ -58,6 +58,57 @@ 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...) +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 pFqweniger(α, β, z; kwds...)::T + end +end +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...)::T + else + return pFqweniger(α, β, z; kwds...)::T + end + elseif p == q + 1 + if ispositivelybalanced(α, β) + #return pFqconformalrationalgeneralizedeuler(α, β, z; kwds...)::T + return pFqconformalrational(α, β, z; kwds...)::T + else + return pFqconformalrational(α, β, z; kwds...)::T + end + else + return pFqweniger(α, β, z; kwds...)::T + end +end + +#= """ Remove parameters common to `α` and `β`, since `(x)ₖ/(x)ₖ ≡ 1` for any shared `x`. Cancelling them avoids loss of accuracy (or a `0/0`, @@ -107,19 +158,71 @@ function pFq(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q} 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 +``` +julia> using HypergeometricFunctions, BenchmarkTools + +julia> α = (3,4,5) +(3, 4, 5) + +julia> β = (2,6,7) +(2, 6, 7) + +julia> z = 0.123 +0.123 + +julia> @time pFq(α, β, z) + 0.253318 seconds (363.21 k allocations: 22.849 MiB, 10.90% gc time, 99.98% compilation time) +1.091834662493663 + +julia> @btime pFq($α, $β, $z) + 1.605 μs (25 allocations: 624 bytes) +1.091834662493663 + +julia> α = (1,1,1,1,1,3,4,5) +(1, 1, 1, 1, 1, 3, 4, 5) + +julia> β = (1,1,1,1,1,2,6,7) +(1, 1, 1, 1, 1, 2, 6, 7) + +julia> @btime pFq($α, $β, $z) + 3.035 μs (50 allocations: 1.22 KiB) +1.091834662493663 +``` + +```julia +julia> using HypergeometricFunctions, BenchmarkTools + +julia> α = (3,4,5) +(3, 4, 5) + +julia> β = (2,6,7) +(2, 6, 7) + +julia> z = 0.123 +0.123 + +julia> @time pFq(α, β, z) + 0.463688 seconds (1.07 M allocations: 68.306 MiB, 5.75% gc time, 99.99% compilation time) +1.091834662493663 + +julia> @btime pFq($α, $β, $z) + 76.297 ns (0 allocations: 0 bytes) +1.091834662493663 + +julia> α = (1,1,1,1,1,3,4,5) +(1, 1, 1, 1, 1, 3, 4, 5) + +julia> β = (1,1,1,1,1,2,6,7) +(1, 1, 1, 1, 1, 2, 6, 7) + +julia> @btime pFq($α, $β, $z) + 317.771 ns (32 allocations: 1.19 KiB) +1.091834662493663 +``` +=# + +pFq(α::AbstractVector, β::AbstractVector, z; kwds...) = pFq(Tuple(α), Tuple(β), z; kwds...) """ Compute the generalized hypergeometric function `₃F₂(a₁, 1, 1, b₁, 2, z)`. diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl index e4ccc32..961c218 100644 --- a/src/specialfunctions.jl +++ b/src/specialfunctions.jl @@ -82,12 +82,12 @@ const ρϵ = 0.71 struct ℕ end Base.in(n::Integer, ::Type{ℕ}) = n > 0 -Base.in(n::Number, ::Type{ℕ}) = isinteger(n) && round(Int, n) ∈ ℕ +Base.in(n::Number, ::Type{ℕ}) = isinteger(n) && Int(real(n)) ∈ ℕ struct ℕ₀ end Base.in(n::Integer, ::Type{ℕ₀}) = n ≥ 0 -Base.in(n::Number, ::Type{ℕ₀}) = isinteger(n) && round(Int, n) ∈ ℕ₀ +Base.in(n::Number, ::Type{ℕ₀}) = isinteger(n) && Int(real(n)) ∈ ℕ₀ struct ℤ end @@ -96,8 +96,56 @@ Base.in(n::Number, ::Type{ℤ}) = isinteger(n) abeqcd(a, b, cd) = isequal(a, b) && isequal(b, cd) abeqcd(a, b, c, d) = isequal(a, c) && isequal(b, d) -iswellpoised(a, b, c) = real(c - a - b) > 0 -isalmostwellpoised(a, b, c) = real(c - a - b) ≥ 0 +ispositivelybalanced(a, b, c) = real(c - a - b) > 0 +isalmostpositivelybalanced(a, b, c) = real(c - a - b) ≥ 0 + +ispositivelybalanced(α::NTuple{p, Any}, β::NTuple{q, Any}) where {p, q} = real(sum(β)-sum(α)) > 0 +isalmostpositivelybalanced(α::NTuple{p, Any}, β::NTuple{q, Any}) where {p, q} = real(sum(β)-sum(α)) ≥ 0 + +function isacommonparameter(α::NTuple{p, S}, β::NTuple{q, T}) where {p, q, S, T} + for i in 1:p + for j in 1:q + if α[i] == β[j] return true end + end + end + return false +end + +function findacommonparameter(α::NTuple{p, S}, β::NTuple{q, T}) where {p, q, S, T} + for i in 1:p + for j in 1:q + if α[i] == β[j] return (i, j) end + end + end + return (p+1, q+1) +end + +""" +Cancel a common parameter in `α` 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. +""" +function cancelacommonparameter(α::NTuple{p, S}, β::NTuple{q, T}) where {p, q, S, T} + i, j = findacommonparameter(α, β) + γ = ntuple(k -> k < i ? α[k] : α[k+1], Val(p-1)) + δ = ntuple(k -> k < j ? β[k] : β[k+1], Val(q-1)) + return γ, δ +end + +function isterminating(α::NTuple{p, Any}) where p + for i in 1:p + if -α[i] ∈ ℕ₀ return true end + end + return false +end + +function getterminatinglength(α::NTuple{p, Any}) where p + kmax = KMAX + for i in 1:p + if -α[i] ∈ ℕ₀ kmax = min(kmax, Int(-α[i])) end + end + return kmax +end absarg(z) = abs(angle(z)) @@ -391,10 +439,9 @@ function AInf(a, b, c, w, m::Int, ϵ) ret end -function BInf(a, b, c, win, m::Int, ϵ) - T = promote_type(typeof(a), typeof(b), typeof(c), typeof(win), typeof(m), typeof(ϵ)) - w = win - βₙ, γₙ = recInfβ₀(a, b, c, win, m, ϵ)*one(w), recInfγ₀(a, b, c, win, m, ϵ)*w +function BInf(a, b, c, w, m::Int, ϵ) + T = promote_type(typeof(a), typeof(b), typeof(c), typeof(w), typeof(m), typeof(ϵ)) + βₙ, γₙ = recInfβ₀(a, b, c, w, m, ϵ)*one(w), recInfγ₀(a, b, c, w, m, ϵ)*w ret, n = βₙ, 0 while abs(βₙ) > 8abs(ret)*eps(real(T)) || 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-ϵ)) @@ -406,6 +453,9 @@ function BInf(a, b, c, win, m::Int, ϵ) end function _₂F₁Inf(a, b, c, z) + if real(b-a) < 0 + return _₂F₁Inf(b, a, c, z) + end m = round(Int, real(b-a)) ϵ = b-a-m w = inv(z) @@ -540,6 +590,36 @@ function pFqmaclaurin(a::NTuple{p, S}, b::NTuple{q, U}, z::V; kmax::Int = KMAX) return S₁ end +function pFqmaclaurin(α::NTuple{p, Any}, β::NTuple{q, Any}, z, kmax; kwds...) where {p, q} + T1 = isempty(α) ? Any : mapreduce(typeof, promote_type, α) + T2 = isempty(β) ? Any : mapreduce(typeof, promote_type, β) + pFqmaclaurin(T1.(α), T2.(β), z, kmax; kwds...) +end + +function pFqmaclaurin(a::NTuple{p, S}, b::NTuple{q, U}, z::V, kmax) where {p, q, S, U, V} + T = promote_type(eltype(a), eltype(b), V) + S₀, S₁, k = one(T), one(T)+prod(a)*z/prod(b), 1 + for k in 1:kmax-1 + rₖ = z/(k+one(T)) + for i=1:p rₖ *= a[i]+k end + for i=1:q rₖ /= b[i]+k end + S₀, S₁ = S₁, S₁+(S₁-S₀)*rₖ + end + return S₁ +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 + function continuedfraction(v::V, u::U, rtol::T) where {V<:Function, U<:Function, T<:Number} n = 4 a, b = continuedfraction(v, u, n, 2n) diff --git a/test/runtests.jl b/test/runtests.jl index e4875a0..e19280c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,11 @@ using HypergeometricFunctions, Gamma, Test import LinearAlgebra: norm -import HypergeometricFunctions: iswellpoised, isalmostwellpoised, M, U, - pochhammer, unsafe_gamma, _₂F₁general, _₂F₁general2, - pFqdrummond, pFqweniger, pFq2string, +import HypergeometricFunctions: M, U, pochhammer, unsafe_gamma, + _₂F₁general, _₂F₁general2, pFqdrummond, pFqweniger, p, r, pFqconformalpolynomial, pFqconformalrational, - cancelcommonparameters + ispositivelybalanced, isalmostpositivelybalanced, + isacommonparameter, cancelacommonparameter, + isterminating, getterminatinglength, pFq2string const rtol = 1.0e-3 const NumberType = Float64 @@ -100,7 +101,7 @@ end for (a, b, c, z) in ((1, 2, 3, 1), (1, 2, 4, 1), (3, 5, 7, 1.0), (1, 8537, 6042, 1.0)) twoFone = _₂F₁(a, b, c, z) twoFonep = _₂F₁(a, b, c, z; method = :positive) - if iswellpoised(a, b, c) + if ispositivelybalanced(a, b, c) @test isfinite(twoFone) @test isfinite(twoFonep) else @@ -109,7 +110,7 @@ end end end for (a, b, c, z) in ((1, 2+im, 3, 1), (1, 2+im, 3+im, 1), (1, 2 + 1im, 2, 1)) - if isalmostwellpoised(a, b, c) + if isalmostpositivelybalanced(a, b, c) if a+b==c @test isinf(_₂F₁(a, b, c, z)) else @@ -632,29 +633,31 @@ end ((1e4, 2e4, 3.0), (4.0, 5.0), 0.00001), ((1.1, 1.2, 1.3), (2.6, 2.7), 0.5), ((1.1, 1.2, 1.3, 1.4), (2.6, 2.7, 2.8), 0.5), - ((1.1, 1.2, 1.3, 1.4, 1.5), (2.6, 2.7, 2.8, 2.9), 2.0+im)) + ((1.1, 1.2, 1.3, 1.4, 1.5), (2.6, 2.7, 2.8, 2.9), 2.0+im), + ((1.9, 2.25, 4.5), (1.17, 2.29), 1.0+(1/8)*im), + ((1.9, 2.25, -4.5), (1.17, 2.29), 1.0+(1/8)*im)) @test pFqconformalrational(α, β, z, Val(2)) ≈ pFqconformalrational(α, β, z, Val(3)) ≈ pFqconformalrational(α, β, z, Val(4)) ≈ pFq(α, β, z) end end end -@testset "Cancellation of common parameters (issue #99)" begin +@testset "Canceling common parameters (issue #99)" begin # pFq((2,3,5,7,11), (2,3,5,7), 2) reduces to ₁F₀(11;;2) = (1-2)^-11 = -1, - @test cancelcommonparameters((2.0, 3.0, 5.0, 7.0, 11.0), (2.0, 3.0, 5.0, 7.0)) == ((11.0,), ()) - @test cancelcommonparameters((2, 3, 5, 7, 11), (2, 3, 5, 7)) == ((11,), ()) - @test cancelcommonparameters((1.0, 2.0), (3.0, 4.0)) == ((1.0, 2.0), (3.0, 4.0)) # no common parameters - @test cancelcommonparameters((1.0, 2.0, 2.0), (2.0,)) == ((1.0, 2.0), ()) # only one instance of a repeated value cancels - - # direct tests of cancelcommonparameters with known input/output pairs - @test cancelcommonparameters((), ()) == ((), ()) - @test cancelcommonparameters((1, 2), ()) == ((1, 2), ()) # empty β - @test cancelcommonparameters((), (1, 2)) == ((), (1, 2)) # empty α - @test cancelcommonparameters((2, 2, 3), (2, 2)) == ((3,), ()) # both repeats in β cancel - @test cancelcommonparameters((1, 2, 3), (3, 2, 1)) == ((), ()) # order-independent, fully cancels - @test cancelcommonparameters((2, 3.0), (2.0,)) == ((3.0,), ()) # cross-type equality (Int 2 == Float64 2.0) still cancels - @test cancelcommonparameters((1, 2, 2, 2), (2, 2)) == ((1, 2), ()) # only as many repeats as present in β cancel - @test cancelcommonparameters((-4, 3, 5), (-4, 9)) == ((3, 5), (9,)) # negative-integer shared value cancels too - @test cancelcommonparameters((2+3im, 5, 7), (2+3im, 11)) == ((5, 7), (11,)) # complex shared value cancels too + α = (2.0, 3.0, 5.0, 7.0, 11.0) + β = (2, 3, 5, 7) + @test isacommonparameter(α, β) + α, β = cancelacommonparameter(α, β) + @test (α, β) == ((3.0, 5.0, 7.0, 11.0), (3, 5, 7)) + @test isacommonparameter(α, β) + α, β = cancelacommonparameter(α, β) + @test (α, β) == ((5.0, 7.0, 11.0), (5, 7)) + @test isacommonparameter(α, β) + α, β = cancelacommonparameter(α, β) + @test (α, β) == ((7.0, 11.0), (7, )) + @test isacommonparameter(α, β) + α, β = cancelacommonparameter(α, β) + @test (α, β) == ((11.0, ), ()) + @test !isacommonparameter(α, β) @test pFq((2, 3, 5, 7, 11), (2, 3, 5, 7), 2) == -1 @test pFq((2.0, 3.0, 5.0, 7.0, 11.0), (2.0, 3.0, 5.0, 7.0), 2) ≈ -1 @@ -678,7 +681,7 @@ end # (4,3) -> (3,2), residual call still dispatches to the generic pFq method @test pFq((1.5, 2.5, 3.5, 9.5), (1.5, 4.5, 6.5), 0.3) ≈ 2.8474959586954557115949150507520688701294357489327 rtol=1e-14 - @test pFq((1.5, 2.5, 3.5, 9.5), (1.5, 4.5, 6.5), 2.0+0.3im) ≈ 0.36311499245824637618815668424357819249646478168082 + 0.044468831058723770609686622207959873127727841990114im rtol=1e-13 + @test pFq((1.5, 2.5, 3.5, 9.5), (1.5, 4.5, 6.5), 2.0+0.3im) ≈ 0.36311499245824637618815668424357819249646478168082 + 0.044468831058723770609686622207959873127727841990114im rtol=1e-10 @test setprecision(BigFloat, 256) do pFq((big"1.5", big"2.5", big"3.5", big"9.5"), (big"1.5", big"4.5", big"6.5"), big"0.3") end ≈ big"2.8474959586954557115949150507520688701294357489327" rtol=big"1e-45" @@ -695,6 +698,19 @@ end @test pFq((2+3im, 5.0, 7.0), (2+3im, 11.0), 2.0+0.4im) ≈ pFq((5.0, 7.0), (11.0,), 2.0+0.4im) end +@testset "Terminating pFq" begin + @test isterminating((-1,2,3)) + @test isterminating((1,-2.0,3)) + @test isterminating((1,2//1,-3.f0+0im)) + @test getterminatinglength((-1,2,3)) == 1 + @test getterminatinglength((-1,-2,3)) == 1 + @test getterminatinglength((-1,-2,-3)) == 1 + @test getterminatinglength((1,-2,-3)) == 2 + @test getterminatinglength((1,2,-3)) == 3 + @test getterminatinglength((-1,2,-3)) == 1 + @test getterminatinglength((1,-2,3)) == 2 +end + @testset "M" begin @test M(-3, -3, 0.5) ≡ exp(0.5) @test M(0, -1, 10) ≡ 1.0 From 6d00e0cbf97b5ea91751e380d79efb11e155fb3d Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Wed, 29 Jul 2026 11:45:07 -0500 Subject: [PATCH 2/3] delete commented code --- src/generalized.jl | 114 --------------------------------------------- 1 file changed, 114 deletions(-) diff --git a/src/generalized.jl b/src/generalized.jl index 70aba3b..5f40899 100644 --- a/src/generalized.jl +++ b/src/generalized.jl @@ -108,120 +108,6 @@ function pFq(α::NTuple{p, T1}, β::NTuple{q, T2}, z; kwds...) where {p, q, T1, end end -#= -""" -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′...), β″ - else - return cancelcommonparameters(resta, β′) - 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...) - end - z = float(z) - if p ≤ q - if real(z) > 0 - return pFqmaclaurin(α, β, z; kwds...) - else - return pFqweniger(α, β, z; kwds...) - end - elseif p == q + 1 - if abs(z) ≤ ρ - return pFqmaclaurin(α, β, z; kwds...) - else - return pFqweniger(α, β, z; kwds...) - end - else - return pFqweniger(α, β, z; kwds...) - end -end - - -``` -julia> using HypergeometricFunctions, BenchmarkTools - -julia> α = (3,4,5) -(3, 4, 5) - -julia> β = (2,6,7) -(2, 6, 7) - -julia> z = 0.123 -0.123 - -julia> @time pFq(α, β, z) - 0.253318 seconds (363.21 k allocations: 22.849 MiB, 10.90% gc time, 99.98% compilation time) -1.091834662493663 - -julia> @btime pFq($α, $β, $z) - 1.605 μs (25 allocations: 624 bytes) -1.091834662493663 - -julia> α = (1,1,1,1,1,3,4,5) -(1, 1, 1, 1, 1, 3, 4, 5) - -julia> β = (1,1,1,1,1,2,6,7) -(1, 1, 1, 1, 1, 2, 6, 7) - -julia> @btime pFq($α, $β, $z) - 3.035 μs (50 allocations: 1.22 KiB) -1.091834662493663 -``` - -```julia -julia> using HypergeometricFunctions, BenchmarkTools - -julia> α = (3,4,5) -(3, 4, 5) - -julia> β = (2,6,7) -(2, 6, 7) - -julia> z = 0.123 -0.123 - -julia> @time pFq(α, β, z) - 0.463688 seconds (1.07 M allocations: 68.306 MiB, 5.75% gc time, 99.99% compilation time) -1.091834662493663 - -julia> @btime pFq($α, $β, $z) - 76.297 ns (0 allocations: 0 bytes) -1.091834662493663 - -julia> α = (1,1,1,1,1,3,4,5) -(1, 1, 1, 1, 1, 3, 4, 5) - -julia> β = (1,1,1,1,1,2,6,7) -(1, 1, 1, 1, 1, 2, 6, 7) - -julia> @btime pFq($α, $β, $z) - 317.771 ns (32 allocations: 1.19 KiB) -1.091834662493663 -``` -=# - pFq(α::AbstractVector, β::AbstractVector, z; kwds...) = pFq(Tuple(α), Tuple(β), z; kwds...) """ From 6a765502656b2a9f8eb769f14df5fbb90409e7e9 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Wed, 29 Jul 2026 12:06:00 -0500 Subject: [PATCH 3/3] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 354589c..bfc57ba 100644 --- a/Project.toml +++ b/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"