Skip to content

Commit

Permalink
add some methods that are ratios of rationals
Browse files Browse the repository at this point in the history
reorganize where float(z) is called

change $kmax to string(kmax)

KMAX = 2^20 instead of 10^6
  • Loading branch information
MikaelSlevinsky committed Jul 6, 2023
1 parent 089e208 commit df46b74
Show file tree
Hide file tree
Showing 9 changed files with 216 additions and 203 deletions.
4 changes: 3 additions & 1 deletion LICENSE
Original file line number Diff line number Diff line change
@@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/HypergeometricFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
1 change: 1 addition & 0 deletions src/confluent.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
131 changes: 54 additions & 77 deletions src/drummond.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand Down
7 changes: 4 additions & 3 deletions src/gauss.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 ℕ₀
Expand Down Expand Up @@ -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 ℕ₀
Expand Down
15 changes: 8 additions & 7 deletions src/generalized.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading

0 comments on commit df46b74

Please sign in to comment.