diff --git a/LICENSE b/LICENSE index 8d19d47..999b425 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,8 @@ MIT License -Copyright (c) 2018 JuliaApproximation +Copyright (c) 2018-2023 Richard Mikael Slevinsky and other contributors: + +https://github.com/JuliaMath/HypergeometricFunctions.jl/graphs/contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/src/HypergeometricFunctions.jl b/src/HypergeometricFunctions.jl index b3965be..9866abc 100644 --- a/src/HypergeometricFunctions.jl +++ b/src/HypergeometricFunctions.jl @@ -9,7 +9,7 @@ using DualNumbers, LinearAlgebra, SpecialFunctions export _₁F₁, _₂F₁, _₃F₂, pFq -const KMAX = 1_000_000 +const KMAX = 1_048_576 include("specialfunctions.jl") include("gauss.jl") diff --git a/src/confluent.jl b/src/confluent.jl index 6cbc0d7..40a1497 100644 --- a/src/confluent.jl +++ b/src/confluent.jl @@ -4,6 +4,7 @@ Compute Kummer's confluent hypergeometric function `M(a, b, z) = ₁F₁(a; b; z)`. """ function _₁F₁(a, b, z; kwds...) + z = float(z) if isequal(a, b) # 13.6.1 return exp(z) elseif -b ∈ ℕ diff --git a/src/drummond.jl b/src/drummond.jl index 6d3e605..7403553 100644 --- a/src/drummond.jl +++ b/src/drummond.jl @@ -6,25 +6,19 @@ function pFqdrummond(::Tuple{}, ::Tuple{}, z::T; kmax::Int = KMAX) where T if norm(z) < eps(real(T)) return one(T) end - ζ = inv(z) - Nlo = ζ - Dlo = ζ - Tlo = Nlo/Dlo - Nhi = (2ζ - 1)*Nlo + 2ζ - Dhi = (2ζ - 1)*Dlo - Thi = Nhi/Dhi - k = 1 - Nhi, Nlo = ((k+2)*ζ-1)*Nhi + k*ζ*Nlo + ζ, Nhi - Dhi, Dlo = ((k+2)*ζ-1)*Dhi + k*ζ*Dlo, Dhi - Thi, Tlo = Nhi/Dhi, Thi - k += 1 + μlo = one(z) + μhi = inv(2-z) + Tlo = one(z) + Thi = Tlo + 2z*μhi + μhi, μlo = inv(3-z + z*μhi), μhi + Thi, Tlo = ((3-z)*Thi + z*(Tlo+z)*μlo)*μhi, Thi + k = 2 while k < kmax && errcheck(Tlo, Thi, 8eps(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 + μhi, μlo = inv(k+2-z + k*z*μhi), μhi + Thi, Tlo = ((k+2-z)*Thi + k*z*Tlo*μlo)*μhi, Thi k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : Tlo end @@ -48,16 +42,15 @@ function pFqdrummond(α::Tuple{T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) where {T end Nhi /= α+1 Dhi /= α+1 - k = 1 - Nhi, Nlo = ((k+2)*ζ-(α+2k+1))*Nhi + k*(ζ-1)*Nlo + ζ, Nhi - Dhi, Dlo = ((k+2)*ζ-(α+2k+1))*Dhi + k*(ζ-1)*Dlo, Dhi + Nhi, Nlo = (3ζ-(α+3))*Nhi + (ζ-1)*Nlo + ζ, Nhi + Dhi, Dlo = (3ζ-(α+3))*Dhi + (ζ-1)*Dlo, Dhi Thi, Tlo = Nhi/Dhi, Thi - if norm(α+k+1) < eps(absα+k+1) + if norm(α+2) < eps(absα+2) return Thi end - Nhi /= α+k+1 - Dhi /= α+k+1 - k += 1 + Nhi /= α+2 + Dhi /= α+2 + k = 2 while k < kmax && errcheck(Tlo, Thi, 8eps(real(T))) Nhi, Nlo = ((k+2)*ζ-(α+2k+1))*Nhi + k*(ζ-1)*Nlo, Nhi Dhi, Dlo = ((k+2)*ζ-(α+2k+1))*Dhi + k*(ζ-1)*Dlo, Dhi @@ -69,7 +62,7 @@ function pFqdrummond(α::Tuple{T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) where {T Dhi /= α+k+1 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : Tlo end @@ -101,7 +94,7 @@ function pFqdrummond(::Tuple{}, β::Tuple{T1}, z::T2; kmax::Int = KMAX) where {T Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{1}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{1}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo end @@ -114,39 +107,32 @@ function pFqdrummond(α::Tuple{T1, T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) wher if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ) return one(T) end - ζ = inv(z) - Nlo = ζ/(α*β) - Dlo = ζ/(α*β) - Tlo = Nlo/Dlo - Nmid = (2ζ-(α+1)*(β+1))*Nlo + 2ζ - Dmid = (2ζ-(α+1)*(β+1))*Dlo - Tmid = Nmid/Dmid + μlo = T(α*β) + Tlo = one(T) + μmid = inv(2-(α+1)*(β+1)*z) + Tmid = Tlo + 2z*μmid*μlo if norm((α+1)*(β+1)) < eps((absα+1)*(absβ+1)) return Tmid end - Nmid /= (α+1)*(β+1) - Dmid /= (α+1)*(β+1) - Nhi = (3ζ-(α+2)*(β+2)-(α+β+3))*Nmid - (α+β+3-ζ)*Nlo + ζ - Dhi = (3ζ-(α+2)*(β+2)-(α+β+3))*Dmid - (α+β+3-ζ)*Dlo - Thi = Nhi/Dhi + μmid *= (α+1)*(β+1) + μhi = inv((3-((α+2)*(β+2)+(α+β+3))*z) + (1-(α+β+3)*z)*z*μmid) + Thi = ((3-((α+2)*(β+2)+(α+β+3))*z)*Tmid + ((1-(α+β+3)*z)*Tlo + z*μlo)*z*μmid)*μhi if norm((α+2)*(β+2)) < eps((absα+2)*(absβ+2)) return Thi end - Nhi /= (α+2)*(β+2) - Dhi /= (α+2)*(β+2) + μhi *= (α+2)*(β+2) k = 2 + z2 = z*z while k < 3 || (k < kmax && errcheck(Tmid, Thi, 8eps(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 + μhi, μmid, μlo = inv((k+2)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z + k*((1-(α+β+3k)*z) - (k-1)*z2*μmid)*z*μhi), μhi, μmid + Thi, Tmid, Tlo = (((k+2)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z)*Thi + k*((1-(α+β+3k)*z)*Tmid - (k-1)*z2*Tlo*μlo)*z*μmid)*μhi, Thi, Tmid if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1)) return Thi end - Nhi /= (α+k+1)*(β+k+1) - Dhi /= (α+k+1)*(β+k+1) + μhi *= (α+k+1)*(β+k+1) k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo end @@ -200,7 +186,7 @@ function pFqdrummond(α::Tuple{T1}, β::Tuple{T2}, z::T3; kmax::Int = KMAX) wher Dhi /= α+k+1 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{1}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{1}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo end @@ -238,7 +224,7 @@ function pFqdrummond(::Tuple{}, β::Tuple{T1, T1}, z::T2; kmax::Int = KMAX) wher Thi, Tmid1, Tmid2, Tlo = Nhi/Dhi, Thi, Tmid1, Tmid2 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{2}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{2}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid1) ? Tmid1 : isfinite(Tmid2) ? Tmid2 : Tlo end @@ -252,59 +238,50 @@ function pFqdrummond(α::Tuple{T1, T1}, β::Tuple{T2}, z::T3; kmax::Int = KMAX) if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ) return one(T) end - ζ = inv(z) - Nlo = γ*ζ/(α*β) - Dlo = γ*ζ/(α*β) - Tlo = Nlo/Dlo - Nmid = ((2)*(γ+1)*ζ-(α+1)*(β+1))*Nlo + 2*(γ+1)*ζ - Dmid = ((2)*(γ+1)*ζ-(α+1)*(β+1))*Dlo - Tmid = Nmid/Dmid + μlo = T(α*β)/T(γ) + Tlo = one(T) + μmid = inv((2)*(γ+1)-(α+1)*(β+1)*z) + Tmid = Tlo + (γ+1)*2z*μmid*μlo if norm((α+1)*(β+1)) < eps((absα+1)*(absβ+1)) return Tmid end - Nmid /= (α+1)*(β+1) - Dmid /= (α+1)*(β+1) - Nhi = ((3)*(γ+2)*ζ-(α+2)*(β+2)-(α+β+3))*Nmid + ((γ+4)*ζ-(α+β+3))*Nlo + (γ+4)*ζ - Dhi = ((3)*(γ+2)*ζ-(α+2)*(β+2)-(α+β+3))*Dmid + ((γ+4)*ζ-(α+β+3))*Dlo - Thi = Nhi/Dhi + μmid *= (α+1)*(β+1) + μhi = inv((3)*(γ+2)-((α+2)*(β+2)+(α+β+3))*z + ((γ+4)-(α+β+3)*z)*z*μmid) + Thi = (((3)*(γ+2)-((α+2)*(β+2)+(α+β+3))*z)*Tmid + (((γ+4)-(α+β+3)*z)*z*Tlo + (γ+4)*z*z*μlo)*μmid)*μhi if norm((α+2)*(β+2)) < eps((absα+2)*(absβ+2)) return Thi end - Nhi /= (α+2)*(β+2) - Dhi /= (α+2)*(β+2) + μhi *= (α+2)*(β+2) k = 2 - 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 + 2ζ, 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 + μlo2 = μlo + μhi, μmid, μlo = inv((k+2)*(γ+k+1)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z + k*(((γ+2k+2)-(α+β+3k)*z)*z + (k-1)*(1-z)*z*z*μmid)*μhi), μhi, μmid + Thi, Tmid, Tlo = (((k+2)*(γ+k+1)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z)*Thi + (k*((γ+2k+2)-(α+β+3k)*z)*z*Tmid + (k*(k-1)*(1-z)*z*z*Tlo + 2z*z*z*μlo2)*μlo)*μmid)*μhi, Thi, Tmid if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1)) return Thi end - Nhi /= (α+k+1)*(β+k+1) - Dhi /= (α+k+1)*(β+k+1) + μhi *= (α+k+1)*(β+k+1) k += 1 while k < kmax && errcheck(Tmid, Thi, 8eps(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 + μhi, μmid, μlo = inv((k+2)*(γ+k+1)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z + k*(((γ+2k+2)-(α+β+3k)*z)*z + (k-1)*(1-z)*z*z*μmid)*μhi), μhi, μmid + Thi, Tmid, Tlo = (((k+2)*(γ+k+1)-((α+k+1)*(β+k+1)+k*(α+β+2k+1))*z)*Thi + k*(((γ+2k+2)-(α+β+3k)*z)*z*Tmid + (k-1)*(1-z)*z*z*Tlo*μlo)*μmid)*μhi, Thi, Tmid if norm((α+k+1)*(β+k+1)) < eps((absα+k+1)*(absβ+k+1)) return Thi end - Nhi /= (α+k+1)*(β+k+1) - Dhi /= (α+k+1)*(β+k+1) + μhi *= (α+k+1)*(β+k+1) k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{1}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{1}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo end # ₘFₙ(α;β;z) -function pFqdrummond(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3; kwds...) where {T1, T2, T3} - pFqdrummond(Tuple(α), Tuple(β), z; kwds...) +function pFqdrummond(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3, args...; kwds...) where {T1, T2, T3} + pFqdrummond(Tuple(α), Tuple(β), z, args...; kwds...) end -function pFqdrummond(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q} +function pFqdrummond(α::NTuple{p, Any}, β::NTuple{q, Any}, z, args...; kwds...) where {p, q} T1 = isempty(α) ? Any : mapreduce(typeof, promote_type, α) T2 = isempty(β) ? Any : mapreduce(typeof, promote_type, β) - pFqdrummond(T1.(α), T2.(β), z; kwds...) + pFqdrummond(T1.(α), T2.(β), z, args...; kwds...) end function pFqdrummond(α::NTuple{p, T1}, β::NTuple{q, T2}, z::T3; kmax::Int = KMAX) where {p, q, T1, T2, T3} T = promote_type(eltype(α), eltype(β), T3) @@ -384,7 +361,7 @@ function pFqdrummond(α::NTuple{p, T1}, β::NTuple{q, T2}, z::T3; kmax::Int = KM end Q[q+2] = t end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(R[r+2]) ? R[r+2] : R[r+1] end diff --git a/src/gauss.jl b/src/gauss.jl index c59be23..6c96e8b 100644 --- a/src/gauss.jl +++ b/src/gauss.jl @@ -4,6 +4,7 @@ Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)`. """ function _₂F₁(a, b, c, z; method::Symbol = :general, kwds...) + z = float(z) if real(b) < real(a) return _₂F₁(b, a, c, z; method = method, kwds...) # ensure a ≤ b elseif isequal(a, c) # 1. 15.4.6 @@ -60,7 +61,7 @@ Compute the Gauss hypergeometric function `₂F₁(a, b; c; z)` with positive pa function _₂F₁positive(a, b, c, z; kwds...) @assert a > 0 && b > 0 && c > 0 && 0 ≤ z ≤ 1 if z == 1 - return _₂F₁argument_unity(a, b, c, float(z); kwds...) + return _₂F₁argument_unity(a, b, c, z; kwds...) else return _₂F₁maclaurin(a, b, c, z; kwds...) end @@ -74,7 +75,7 @@ N. Michel and M. V. Stoitsov, Fast computation of the Gauss hypergeometric funct """ function _₂F₁general(a, b, c, z; kwds...) if z == 1 - return _₂F₁argument_unity(a, b, c, float(z); kwds...) + return _₂F₁argument_unity(a, b, c, z; kwds...) elseif real(b) < real(a) return _₂F₁general(b, a, c, z; kwds...) elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀ @@ -105,7 +106,7 @@ J. W. Pearson, S. Olver and M. A. Porter, Numerical methos for the computation o function _₂F₁general2(a, b, c, z) T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z)) if z == 1 - return _₂F₁argument_unity(a, b, c, float(z)) + return _₂F₁argument_unity(a, b, c, z) elseif real(b) < real(a) return _₂F₁general2(b, a, c, z) elseif abs(z) ≤ ρ || -a ∈ ℕ₀ || -b ∈ ℕ₀ diff --git a/src/generalized.jl b/src/generalized.jl index 9541d12..cf550ae 100644 --- a/src/generalized.jl +++ b/src/generalized.jl @@ -9,24 +9,25 @@ pFq pFq(::Tuple{}, ::Tuple{}, z; kwds...) = exp(z) pFq(α::NTuple{1}, ::Tuple{}, z; kwds...) = exp(-α[1]*log1p(-z)) -pFq(α::NTuple{1}, β::NTuple{1}, z; kwds...) = _₁F₁(α[1], β[1], float(z); kwds...) -pFq(α::NTuple{2, Any}, β::NTuple{1}, z; kwds...) = _₂F₁(α[1], α[2], β[1], float(z); kwds...) +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} + z = float(z) if p ≤ q if real(z) ≥ 0 - return pFqmaclaurin(α, β, float(z); kwds...) + return pFqmaclaurin(α, β, z; kwds...) else - return pFqweniger(α, β, float(z); kwds...) + return pFqweniger(α, β, z; kwds...) end elseif p == q + 1 if abs(z) ≤ ρ - return pFqmaclaurin(α, β, float(z); kwds...) + return pFqmaclaurin(α, β, z; kwds...) else - return pFqweniger(α, β, float(z); kwds...) + return pFqweniger(α, β, z; kwds...) end else - return pFqweniger(α, β, float(z); kwds...) + return pFqweniger(α, β, z; kwds...) end end diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl index c551075..4585576 100644 --- a/src/specialfunctions.jl +++ b/src/specialfunctions.jl @@ -414,7 +414,7 @@ end function _₂F₁maclaurin(a::Number, b::Number, c::Number, z::Number) - T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z))) + 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₁, 8eps(real(T))) rⱼ = (a+j)*z/(j+1)*(b+j)/(c+j) @@ -425,7 +425,7 @@ function _₂F₁maclaurin(a::Number, b::Number, c::Number, z::Number) end function _₂F₁maclaurinalt(a::Number, b::Number, c::Number, z::Number) - T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z))) + T = promote_type(typeof(a), typeof(b), typeof(c), typeof(z)) C, S, j = one(T), one(T), 0 while abs(C) > 8abs(S)*eps(real(T)) || j ≤ 1 C *= (a+j)/(j+1)*(b+j)/(c+j)*z @@ -436,7 +436,7 @@ function _₂F₁maclaurinalt(a::Number, b::Number, c::Number, z::Number) end function _₂F₁continuation(s::Number, t::Number, c::Number, z₀::Number, z::Number) - T = float(promote_type(typeof(s), typeof(t), typeof(c), typeof(z₀), typeof(z))) + 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₁, 8eps(real(T))) || j ≤ 2 @@ -448,7 +448,7 @@ function _₂F₁continuation(s::Number, t::Number, c::Number, z₀::Number, z:: end function _₂F₁continuationalt(a::Number, c::Number, z₀::Number, z::Number) - T = float(promote_type(typeof(a), typeof(c), typeof(z₀), typeof(z))) + 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₀ @@ -469,7 +469,7 @@ function _₂F₁continuationalt(a::Number, c::Number, z₀::Number, z::Number) end function _₂F₁logsum(a::Number, b::Number, z::Number, w::Number, s::Int) - T = float(promote_type(typeof(a), typeof(b), typeof(z), typeof(w))) + 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) > 8abs(S)*eps(real(T)) || j ≤ 1 @@ -482,7 +482,7 @@ function _₂F₁logsum(a::Number, b::Number, z::Number, w::Number, s::Int) end function _₂F₁logsumalt(a::Number, b::Number, z::Number, w::Number) - T = float(promote_type(typeof(a), typeof(b), typeof(z), typeof(w))) + 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) > 8abs(S)*eps(real(T)) || j ≤ 1 @@ -495,7 +495,7 @@ function _₂F₁logsumalt(a::Number, b::Number, z::Number, w::Number) end function _₂F₁taylor(a::Number, b::Number, c::Number, z::Number) - T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z))) + 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₀ @@ -510,7 +510,7 @@ function _₂F₁taylor(a::Number, b::Number, c::Number, z::Number) end function _₁F₁maclaurin(a::Number, b::Number, z::Number) - T = float(promote_type(typeof(a), typeof(b), typeof(z))) + T = promote_type(typeof(a), typeof(b), typeof(z)) S₀, S₁, j = one(T), one(T)+a*z/b, 1 while errcheck(S₀, S₁, 8eps(real(T))) rⱼ = (a+j)*z/((b+j)*(j+1)) @@ -527,7 +527,7 @@ function pFqmaclaurin(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where end function pFqmaclaurin(a::NTuple{p, S}, b::NTuple{q, U}, z::V; kmax::Int = KMAX) where {p, q, S, U, V} - T = float(promote_type(eltype(a), eltype(b), V)) + T = promote_type(eltype(a), eltype(b), V) S₀, S₁, k = one(T), one(T)+prod(a)*z/prod(b), 1 while k < kmax && errcheck(S₀, S₁, 8eps(real(T))) rₖ = z/(k+one(T)) @@ -536,7 +536,7 @@ function pFqmaclaurin(a::NTuple{p, S}, b::NTuple{q, U}, z::V; kmax::Int = KMAX) S₀, S₁ = S₁, S₁+(S₁-S₀)*rₖ k += 1 end - k < kmax || @warn "Maclaurin approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum degree of $kmax." + k < kmax || @warn "Maclaurin approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum degree of "*string(kmax)*"." return S₁ end diff --git a/src/weniger.jl b/src/weniger.jl index 25e4ed8..a8b2c7e 100644 --- a/src/weniger.jl +++ b/src/weniger.jl @@ -6,22 +6,19 @@ function pFqweniger(::Tuple{}, ::Tuple{}, z::T; kmax::Int = KMAX) where T if norm(z) < eps(real(T)) return one(T) end - ζ = inv(z) - Nlo = ζ - Dlo = ζ - Tlo = Nlo/Dlo - Nhi = (2ζ - 1)*Nlo + 2ζ - Dhi = (2ζ - 1)*Dlo - Thi = Nhi/Dhi + μlo = one(T) + μhi = inv(2-z) + Rlo = one(T) + Rhi = Rlo + 2z*μhi k = 1 - while k < kmax && errcheck(Tlo, Thi, 8eps(real(T))) - Nhi, Nlo = (4k+2)*ζ*Nhi + Nlo, Nhi - Dhi, Dlo = (4k+2)*ζ*Dhi + Dlo, Dhi - Thi, Tlo = Nhi/Dhi, Thi + z2 = z*z + while k < kmax && errcheck(Rlo, Rhi, 8eps(real(T))) + μhi, μlo = inv(4k+2 + z2*μhi), μhi + Rhi, Rlo = ((4k+2)*Rhi + z2*Rlo*μlo)*μhi, Rhi k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." - return isfinite(Thi) ? Thi : Tlo + k < kmax || @warn "Rational approximation to "*pFq2string(Val{0}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." + return isfinite(Rhi) ? Rhi : Rlo end # ₁F₀(α;z), γ = 2. @@ -32,32 +29,27 @@ function pFqweniger(α::Tuple{T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) where {T1 if norm(z) < eps(real(T)) || norm(α) < eps(absα) return one(T) end - ζ = inv(z) - Nlo = ζ/α - Dlo = ζ/α - Tlo = Nlo/Dlo - Nhi = (2ζ - (α+1))*Nlo + 2ζ - Dhi = (2ζ - (α+1))*Dlo - Thi = Nhi/Dhi + μlo = α + μhi = inv(2-(α+1)*z) + Rlo = one(T) + Rhi = Rlo + 2z*μhi*μlo if norm(α+1) < eps(absα+1) - return Thi + return Rhi end - Nhi /= α+1 - Dhi /= α+1 + μhi *= α+1 k = 1 - while k < kmax && errcheck(Tlo, Thi, 8eps(real(T))) - Nhi, Nlo = (2k+1)*(2ζ-1)*Nhi - (k-α)*Nlo, Nhi - Dhi, Dlo = (2k+1)*(2ζ-1)*Dhi - (k-α)*Dlo, Dhi - Thi, Tlo = Nhi/Dhi, Thi + z2 = z*z + while k < kmax && errcheck(Rlo, Rhi, 8eps(real(T))) + μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi + Rhi, Rlo = ((2k+1)*(2-z)*Rhi - (k-α)*z2*Rlo*μlo)*μhi, Rhi if norm(α+k+1) < eps(absα+k+1) - return Thi + return Rhi end - Nhi /= α+k+1 - Dhi /= α+k+1 + μhi *= α+k+1 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." - return isfinite(Thi) ? Thi : Tlo + k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." + return isfinite(Rhi) ? Rhi : Rlo end # ₂F₀(α,β;z), algorithm γ = 2. @@ -69,54 +61,106 @@ function pFqweniger(α::Tuple{T1, T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) where if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ) return one(T) end + μlo = T(α*β) + Rlo = one(T) + a0 = T((α+1)*(β+1)) + μmid = inv(2-a0*z) + Rmid = Rlo + 2z*μmid*μlo + if norm(a0) < eps((absα+1)*(absβ+1)) + return Rmid + end + μmid *= a0 + k = 1 + a0 = T((α+2)*(β+2)) + t0 = 6-(6-3*α*β)*z + t1 = 6-2*T(2*α*β+α+β-1)*z + μhi = inv(t0 - t1*z*μmid) + Rhi = (t0*Rmid - (t1*Rlo + 6*z*μlo)*z*μmid)*μhi + if norm(a0) < eps((absα+2)*(absβ+2)) + return Rhi + end + μhi *= a0 + k = 2 + z2 = z*z + while k < 3 || (k < kmax && errcheck(Rmid, Rhi, 8eps(real(T)))) + a0 = T((α+k+1)*(β+k+1)) + t0 = (4k+2)-T((k*(α+β+3k)-(α+1)*(β+1)))*T(2k+1)/T(2k-1)*z + t1 = (4k+2)+T(k*(3k-α-β)-(α+1)*(β+1))*z + a2 = T((α+1-k)*(β+1-k))*T(2k+1)/T(2k-1)*z2 + μhi, μmid, μlo = inv(t0 - (t1 + a2*μmid)*z*μhi), μhi, μmid + Rhi, Rmid, Rlo = (t0*Rhi - (t1*Rmid + a2*Rlo*μlo)*z*μmid)*μhi, Rhi, Rmid + if norm(a0) < eps((absα+k+1)*(absβ+k+1)) + return Rhi + end + μhi *= a0 + k += 1 + end + k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." + return isfinite(Rhi) ? Rhi : isfinite(Rmid) ? Rmid : Rlo +end + +# ₁F₁(α,β;z), algorithm γ = 2. +function pFqweniger(α::Tuple{T1}, β::Tuple{T2}, z::T3; kmax::Int = KMAX) where {T1, T2, T3} + α = α[1] + β = β[1] + T = promote_type(T1, T2, T3) + absα = abs(T(α)) + if norm(z) < eps(real(T)) || norm(α) < eps(absα) + return one(T) + end ζ = inv(z) - Nlo = ζ/(α*β) - Dlo = ζ/(α*β) + Nlo = β*ζ/α + Dlo = β*ζ/α Tlo = Nlo/Dlo - a0 = T((α+1)*(β+1)) - Nmid = (2ζ-a0)*Nlo + 2ζ - Dmid = (2ζ-a0)*Dlo + a0 = T(α+1) + b0 = T(2*(β+1)) + Nmid = (b0*ζ-a0)*Nlo + b0*ζ + Dmid = (b0*ζ-a0)*Dlo Tmid = Nmid/Dmid - if norm(a0) < eps((absα+1)*(absβ+1)) + if norm(a0) < eps(absα+1) return Tmid end Nmid /= a0 Dmid /= a0 k = 1 - a0 = T((α+2)*(β+2)) - #a1 = T(2*(2-(2*α*β+α+β+1))) - t0 = 6ζ-6+3*α*β - t1 = 6ζ-2*T(2*α*β+α+β-1) - Nhi = t0*Nmid - t1*Nlo - 6ζ - Dhi = t0*Dmid - t1*Dlo - #Nhi = -a0*Nmid - a1*(Nmid+Nlo) + 6ζ*(Nmid - Nlo) - 6ζ - #Dhi = -a0*Dmid - a1*(Dmid+Dlo) + 6ζ*(Dmid - Dlo) + a0 = T(α+2) + #a1 = -T(4α+2) + b0 = T(6*(β+2)) + b1 = T(-6*β) + t0 = b0*ζ+3α + t1 = b1*ζ+4α+2 + Nhi = t0*Nmid + t1*Nlo + b1*ζ + Dhi = t0*Dmid + t1*Dlo + #Nhi = -a0*Nmid - a1*(Nmid+Nlo) + ζ*(b0*Nmid + b1*Nlo) + b1*ζ + #Dhi = -a0*Dmid - a1*(Dmid+Dlo) + ζ*(b0*Dmid + b1*Dlo) Thi = Nhi/Dhi - if norm(a0) < eps((absα+2)*(absβ+2)) + if norm(a0) < eps(absα+2) return Thi end Nhi /= a0 Dhi /= a0 k = 2 - while k < 3 || (k < kmax && errcheck(Tmid, Thi, 8eps(real(T)))) - a0 = T((α+k+1)*(β+k+1)) - #a1 = T((2*k*k-(2*α*β+α+β+1)))*T(2k)/T(2k-1) - t0 = (4k+2)*ζ-T((k*(α+β+3k)-(α+1)*(β+1)))*T(2k+1)/T(2k-1) - t1 = (4k+2)*ζ+T(k*(3k-α-β)-(α+1)*(β+1)) - a2 = T((α+1-k)*(β+1-k))*T(2k+1)/T(2k-1) - Nhi, Nmid, Nlo = t0*Nhi - t1*Nmid - a2*Nlo, Nhi, Nmid - Dhi, Dmid, Dlo = t0*Dhi - t1*Dmid - a2*Dlo, Dhi, Dmid - #Nhi, Nmid, Nlo = -a0*Nhi - a1*(Nhi+Nmid) - a2*(Nmid+Nlo) + (4k+2)*ζ*(Nhi-Nmid), Nhi, Nmid - #Dhi, Dmid, Dlo = -a0*Dhi - a1*(Dhi+Dmid) - a2*(Dmid+Dlo) + (4k+2)*ζ*(Dhi-Dmid), Dhi, Dmid + while k < 5 || (k < kmax && errcheck(Tmid, Thi, 8eps(real(T)))) + a0 = T(α+k+1) + #a1 = -T(2α+1)*T(2k)/T(2k-1) + a2 = T(α+1-k)*T(2k+1)/T(2k-1) + b0 = T(4k+2)*T(β+k+1) + b1 = T(4k+2)*T(k-β-1) + t0 = b0*ζ+a2 + t1 = b1*ζ+a0 + Nhi, Nmid, Nlo = t0*Nhi + t1*Nmid - a2*Nlo, Nhi, Nmid + Dhi, Dmid, Dlo = t0*Dhi + t1*Dmid - a2*Dlo, Dhi, Dmid + #Nhi, Nmid, Nlo = -a0*Nhi - a1*(Nhi+Nmid) - a2*(Nmid+Nlo) + ζ*(b0*Nhi + b1*Nmid), Nhi, Nmid + #Dhi, Dmid, Dlo = -a0*Dhi - a1*(Dhi+Dmid) - a2*(Dmid+Dlo) + ζ*(b0*Dhi + b1*Dmid), Dhi, Dmid Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid - if norm(a0) < eps((absα+k+1)*(absβ+k+1)) + if norm(a0) < eps(absα+k+1) return Thi end Nhi /= a0 Dhi /= a0 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{0}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{1}(), Val{1}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo end @@ -130,71 +174,58 @@ function pFqweniger(α::Tuple{T1, T1}, β::Tuple{T2}, z::T3; kmax::Int = KMAX) w if norm(z) < eps(real(T)) || norm(α*β) < eps(absα*absβ) return one(T) end - ζ = inv(z) - Nlo = γ*ζ/(α*β) - Dlo = γ*ζ/(α*β) - Tlo = Nlo/Dlo + μlo = T(α*β)/T(γ) + Rlo = one(T) a0 = T((α+1)*(β+1)) b0 = T(2*(γ+1)) - Nmid = (b0*ζ-a0)*Nlo + b0*ζ - Dmid = (b0*ζ-a0)*Dlo - Tmid = Nmid/Dmid + μmid = inv(b0-a0*z) + Rmid = Rlo + b0*z*μmid*μlo if norm(a0) < eps((absα+1)*(absβ+1)) - return Tmid + return Rmid end - Nmid /= a0 - Dmid /= a0 + μmid *= a0 k = 1 a0 = T((α+2)*(β+2)) - #a1 = T(2*(2-(2*α*β+α+β+1))) b0 = T(6*(γ+2)) b1 = T(-6*γ) - t0 = b0*ζ-6+3*α*β - t1 = b1*ζ+2*T(2*α*β+α+β-1) - Nhi = t0*Nmid + t1*Nlo + b1*ζ - Dhi = t0*Dmid + t1*Dlo - #Nhi = -a0*Nmid - a1*(Nmid+Nlo) + ζ*(b0*Nmid + b1*Nlo) + b1*ζ - #Dhi = -a0*Dmid - a1*(Dmid+Dlo) + ζ*(b0*Dmid + b1*Dlo) - Thi = Nhi/Dhi + t0 = b0-(6-3*α*β)*z + t1 = b1+2*T(2*α*β+α+β-1)*z + μhi = inv(t0 + t1*z*μmid) + Rhi = (t0*Rmid + (t1*Rlo + b1*z*μlo)*z*μmid)*μhi if norm(a0) < eps((absα+2)*(absβ+2)) - return Thi + return Rhi end - Nhi /= a0 - Dhi /= a0 + μhi *= a0 k = 2 - while k < 3 || (k < kmax && errcheck(Tmid, Thi, 8eps(real(T)))) + z2 = z*z + while k < 5 || (k < kmax && errcheck(Rmid, Rhi, 8eps(real(T)))) a0 = T((α+k+1)*(β+k+1)) - #a1 = T((2*k*k-(2*α*β+α+β+1)))*T(2k)/T(2k-1) - a2 = T((α+1-k)*(β+1-k))*T(2k+1)/T(2k-1) + a2 = T((α+1-k)*(β+1-k))*T(2k+1)/T(2k-1)*z2 b0 = T(4k+2)*T(γ+k+1) b1 = T(4k+2)*T(k-γ-1) - t0 = b0*ζ-T((k*(α+β+3k)-(α+1)*(β+1)))*T(2k+1)/T(2k-1) - t1 = b1*ζ-T(k*(3k-α-β)-(α+1)*(β+1)) - Nhi, Nmid, Nlo = t0*Nhi + t1*Nmid - a2*Nlo, Nhi, Nmid - Dhi, Dmid, Dlo = t0*Dhi + t1*Dmid - a2*Dlo, Dhi, Dmid - #Nhi, Nmid, Nlo = -a0*Nhi - a1*(Nhi+Nmid) - a2*(Nmid+Nlo) + ζ*(b0*Nhi + b1*Nmid), Nhi, Nmid - #Dhi, Dmid, Dlo = -a0*Dhi - a1*(Dhi+Dmid) - a2*(Dmid+Dlo) + ζ*(b0*Dhi + b1*Dmid), Dhi, Dmid - Thi, Tmid, Tlo = Nhi/Dhi, Thi, Tmid + t0 = b0-T((k*(α+β+3k)-(α+1)*(β+1)))*T(2k+1)/T(2k-1)*z + t1 = b1-T(k*(3k-α-β)-(α+1)*(β+1))*z + μhi, μmid, μlo = inv(t0 + (t1 - a2*μmid)*z*μhi), μhi, μmid + Rhi, Rmid, Rlo = (t0*Rhi + (t1*Rmid - a2*Rlo*μlo)*z*μmid)*μhi, Rhi, Rmid if norm(a0) < eps((absα+k+1)*(absβ+k+1)) - return Thi + return Rhi end - Nhi /= a0 - Dhi /= a0 + μhi *= a0 k += 1 end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{1}())*" reached the maximum type of ($kmax, $kmax)." - return isfinite(Thi) ? Thi : isfinite(Tmid) ? Tmid : Tlo + k < kmax || @warn "Rational approximation to "*pFq2string(Val{2}(), Val{1}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." + return isfinite(Rhi) ? Rhi : isfinite(Rmid) ? Rmid : Rlo end # ₘFₙ(α;β;z) # γ ∉ ℕ -function pFqweniger(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3; kwds...) where {T1, T2, T3} - pFqweniger(Tuple(α), Tuple(β), z; kwds...) +function pFqweniger(α::AbstractVector{T1}, β::AbstractVector{T2}, z::T3, args...; kwds...) where {T1, T2, T3} + pFqweniger(Tuple(α), Tuple(β), z, args...; kwds...) end -function pFqweniger(α::NTuple{p, Any}, β::NTuple{q, Any}, z; kwds...) where {p, q} +function pFqweniger(α::NTuple{p, Any}, β::NTuple{q, Any}, z, args...; kwds...) where {p, q} T1 = isempty(α) ? Any : mapreduce(typeof, promote_type, α) T2 = isempty(β) ? Any : mapreduce(typeof, promote_type, β) - pFqweniger(T1.(α), T2.(β), z; kwds...) + pFqweniger(T1.(α), T2.(β), z, args...; kwds...) end function pFqweniger(α::NTuple{p, T1}, β::NTuple{q, T2}, z::T3; kmax::Int = KMAX, γ::T4 = 2) where {p, q, T1, T2, T3, T4} T = promote_type(eltype(α), eltype(β), T3, T4) @@ -341,6 +372,6 @@ function pFqweniger(α::NTuple{p, T1}, β::NTuple{q, T2}, z::T3; kmax::Int = KMA Q[r+1] = t1 end end - k < kmax || @warn "Rational approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum type of ($kmax, $kmax)." + k < kmax || @warn "Rational approximation to "*pFq2string(Val{p}(), Val{q}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")." return isfinite(R[r+3]) ? R[r+3] : R[r+2] end diff --git a/test/runtests.jl b/test/runtests.jl index 35099aa..9687069 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -379,13 +379,13 @@ end 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 pFqdrummond((α, β), (), z) ≈ S(pFqdrummond((T(α), T(β)), (), T(z))) atol=atol rtol=rtol - @test pFqweniger((α, β), (), z) ≈ S(pFqweniger((T(α), T(β)), (), T(z))) atol=atol rtol=rtol + @test pFqdrummond((α, β), (), z) ≈ S(pFq((T(α), T(β)), (), T(z))) atol=atol rtol=rtol + @test pFqweniger((α, β), (), z) ≈ S(pFq((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 pFqdrummond((α, β), (), z) ≈ CS(pFqdrummond((T(α), T(β)), (), CT(z))) atol=atol rtol=rtol - @test pFqweniger((α, β), (), z) ≈ CS(pFqweniger((T(α), T(β)), (), CT(z))) atol=atol rtol=rtol + @test pFqdrummond((α, β), (), z) ≈ CS(pFq((T(α), T(β)), (), CT(z))) atol=atol rtol=rtol + @test pFqweniger((α, β), (), z) ≈ CS(pFq((T(α), T(β)), (), CT(z))) atol=atol rtol=rtol end end end