Merge branch 'master' into master

Author: TSGut

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.6"
+version = "0.4.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -21,15 +21,15 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 [compat]
 ArrayLayouts = "1"
 BandedMatrices = "0.17"
-ClassicalOrthogonalPolynomials = "0.9, 0.10"
+ClassicalOrthogonalPolynomials = "0.10.1"
 ContinuumArrays = "0.12, 0.13"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.9"
 InfiniteLinearAlgebra = "0.6.19"
 LazyArrays = "1"
 QuasiArrays = "0.9, 0.10"
-SingularIntegrals = "0.0.2, 0.1"
+SingularIntegrals = "0.1"
 SpecialFunctions = "1.0, 2"
 julia = "1.9"
 

src/SemiclassicalOrthogonalPolynomials.jl

@@ -146,10 +146,10 @@ function semiclassical_jacobimatrix(t, a, b, c)
     iszero(c) && return jacobimatrix(P)
     if isone(c)
         return cholesky_jacobimatrix(x->(t-x),P)
-    elseif isone(c÷2)
+    elseif isone(c/2)
         return qr_jacobimatrix(x->(t-x),P)
     elseif isinteger(c) && c ≥ 0 # reduce other integer c cases to hierarchy
-        return SemiclassicalJacobi.(t, a, b, 0:c)[end].X
+        return SemiclassicalJacobi.(t, a, b, 0:Int(c))[end].X
     else # if c is not an integer, use Lanczos for now
         x = axes(P,1)
         return jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), jacobi(b, a, UnitInterval{T}())))
@@ -170,11 +170,11 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
         return Q.X
     end
 
-    if isone(Δa÷2) && iszero(Δb) && iszero(Δc)  # raising by 2
+    if isone(Δa/2) && iszero(Δb) && iszero(Δc)  # raising by 2
         qr_jacobimatrix(Q.X,Q)
-    elseif iszero(Δa) && isone(Δb÷2) && iszero(Δc)
+    elseif iszero(Δa) && isone(Δb/2) && iszero(Δc)
         qr_jacobimatrix(I-Q.X,Q)
-    elseif iszero(Δa) && iszero(Δb) && isone(Δc÷2)
+    elseif iszero(Δa) && iszero(Δb) && isone(Δc/2)
         qr_jacobimatrix(Q.t*I-Q.X,Q)
     elseif isone(Δa) && iszero(Δb) && iszero(Δc)  # raising by 1
         cholesky_jacobimatrix(Q.X,Q)
@@ -188,7 +188,7 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c)
     elseif c > Q.c
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c)
-    # TODO: Implement lowering via QL or via inverting R
+    # TODO: Implement lowering via QL/reverse Cholesky or via inverting R
     # elseif b < Q.b  # iterative lowering by 1
     #    semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c)
     # elseif c < Q.c
@@ -208,6 +208,7 @@ LanczosPolynomial(P::SemiclassicalJacobi{T}) where T =
 
 gives either a mapped `Jacobi` or `LanczosPolynomial` version of `P`.
 """
+# TODO: Use ConvertedOPs for integer special cases?
 toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : LanczosPolynomial(P)
 
 copy(P::SemiclassicalJacobi) = P
@@ -224,8 +225,6 @@ function summary(io::IO, P::SemiclassicalJacobi)
     print(io, "SemiclassicalJacobi with weight x^$a * (1-x)^$b * ($t-x)^$c")
 end
 
-
-
 jacobimatrix(P::SemiclassicalJacobi) = P.X
 
 """
@@ -260,40 +259,51 @@ function semijacobi_ldiv(P::SemiclassicalJacobi, Q)
     (P \ R) * _p0(R̃) * (R̃ \ Q)
 end
 
+# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q in a single step via decomposition.
+function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
+    (Q == P) && return SymTridiagonal(Ones(∞),Zeros(∞))
+    Δa = Q.a-P.a
+    Δb = Q.b-P.b
+    Δc = Q.c-P.c
+    M = Diagonal(Ones(∞))
+    if iseven(Δa) && iseven(Δb) && iseven(Δc)
+        M = qr(P.X^(Δa÷2)*(I-P.X)^(Δb÷2)*(Q.t*I-P.X)^(Δc÷2)).R
+        return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1
+    elseif isone(Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (1,0,0)
+        M = cholesky(P.X).U
+        return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1    
+    elseif iszero(Δa) && isone(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,1,0)
+        M = cholesky(I-P.X).U
+        return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1
+    elseif iszero(Δa) && iszero(Δb) && isone(Δc) # special case (Δa,Δb,Δc) = (0,0,1)
+        M = cholesky(Q.t*I-P.X).U
+        return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1
+    elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc)
+        M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U
+        return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1
+    else
+        error("Implement")
+    end
+end
+
 # returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q.
 function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
     @assert Q.t ≈ P.t
-    # For polynomial modifications, use Cholesky/QR. Use Lanzos otherwise.
-    M = Diagonal(Ones(∞))
-    (Q == P) && return M
+    (Q == P) && return SymTridiagonal(Ones(∞),Zeros(∞))
     Δa = Q.a-P.a
     Δb = Q.b-P.b
     Δc = Q.c-P.c
-    if iseven(Δa) && iseven(Δb) && iseven(Δc)
-        if !iszero(Δa)
-            M = ApplyArray(*,M,P.X^(Δa÷2))
-        end
-        if !iszero(Δb)
-            M = ApplyArray(*,M,(I-P.X)^(Δb÷2))
-        end
-        if !iszero(Δc)
-            M = ApplyArray(*,M,(Q.t*I-P.X)^(Δc÷2))
-        end
-        M = qr(M).R
-        return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match our normalization choice P_0(x) = 1
-    elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc)
-        if !iszero(Δa)
-            M = ApplyArray(*,M,P.X^(Δa))
-        end
-        if !iszero(Δb)
-            M = ApplyArray(*,M,(I-P.X)^(Δb))
-        end
-        if !iszero(Δc)
-            M = ApplyArray(*,M,(Q.t*I-P.X)^(Δc))
+    if isinteger(Δa) && isinteger(Δb) && isinteger(Δc) # (Δa,Δb,Δc) are integers -> use QR/Cholesky iteratively
+        if ((isone(Δa)||isone(Δa/2)) && iszero(Δb) && iszero(Δc)) || (iszero(Δa) && (isone(Δb)||isone(Δb/2)) && iszero(Δc))  || (iszero(Δa) && iszero(Δb) && (isone(Δc)||isone(Δc/2)))
+            M = semijacobi_ldiv_direct(Q, P)
+        elseif Δa > 0  # iterative modification by 1
+            M = ApplyArray(*,semijacobi_ldiv_direct(Q, SemiclassicalJacobi(Q.t, Q.a-1-iseven(Δa), Q.b, Q.c)),semijacobi_ldiv(SemiclassicalJacobi(Q.t, Q.a-1-iseven(Δa), Q.b, Q.c), P))
+        elseif Δb > 0
+            M = ApplyArray(*,semijacobi_ldiv_direct(Q, SemiclassicalJacobi(Q.t, Q.a, Q.b-1-iseven(Δb), Q.c)),semijacobi_ldiv(SemiclassicalJacobi(Q.t, Q.a, Q.b-1-iseven(Δb), Q.c), P))
+        elseif Δc > 0
+            M = ApplyArray(*,semijacobi_ldiv_direct(Q, SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1-iseven(Δc))),semijacobi_ldiv(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1-iseven(Δc)), P))
         end
-        M = cholesky(Symmetric(M)).U
-        return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match our normalization choice P_0(x) = 1
-    else # fallback option for non-integer weight modification
+    else # fallback to Lancos
         R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
         R̃ = toclassical(R)
         return (P \ R) * _p0(R̃) * (R̃ \ Q)
@@ -338,8 +348,14 @@ function \(w_A::WeightedSemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi)
     wA,A = w_A.args
     wB,B = w_B.args
     @assert wA.t == wB.t == A.t == B.t
-
-    if wA.a == wB.a && wA.b == wB.b && wA.c == wB.c
+    Δa = B.a-A.a
+    Δb = B.b-A.b
+    Δc = B.c-A.c
+
+    if (wA.a == A.a) && (wA.b == A.b) && (wA.c == A.c) && (wB.a == B.a) && (wB.b == B.b) && (wB.c == B.c) && isinteger(A.a) && isinteger(A.b) && isinteger(A.c) && isinteger(B.a) && isinteger(B.b) && isinteger(B.c)
+        k = (A \ SemiclassicalJacobiWeight(A.t,Δa,Δb,Δc))[1]
+        return (ApplyArray(*,Diagonal(Fill(k,∞)),(B \ A)))'
+    elseif wA.a == wB.a && wA.b == wB.b && wA.c == wB.c # fallback to Christoffel–Darboux
         A \ B
     elseif wA.a+1 == wB.a && wA.b == wB.b && wA.c == wB.c
         @assert A.a+1 == B.a && A.b == B.b && A.c == B.c
@@ -383,29 +399,30 @@ massmatrix(P::SemiclassicalJacobi) = Diagonal(Fill(sum(orthogonalityweight(P)),
 end
 
 function ldiv(Q::SemiclassicalJacobi, f::AbstractQuasiVector)
-    if iszero(Q.a) && iszero(Q.b) && Q.c == -one(eltype(Q.t))
-        # TODO: due to a stdlib error this won't work with bigfloat as is
-        T = typeof(Q.t)
+    T = typeof(Q.t)
+    if iszero(Q.a) && iszero(Q.b) && isone(-Q.c) # (0,0,-1) special case
         R = legendre(zero(T)..one(T))
         B = neg1c_tolegendre(Q.t)
         return (B \ (R \ f))
+    elseif isinteger(Q.a) && isinteger(Q.b) && isinteger(Q.c) # (a,b,c) are integers -> use QR/Cholesky
+        R̃ = Normalized(jacobi(Q.b, Q.a, UnitInterval{T}()))
+        return (Q \ SemiclassicalJacobi(Q.t, Q.a, Q.b, 0)) *  _p0(R̃) * (R̃ \ f)
+    else # fallback to Lanzcos
+        R̃ = toclassical(SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1)))
+        return (Q \ R̃) * (R̃ \ f)
     end
-    R = SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1))
-    R̃ = toclassical(R)
-    return (Q \ R̃) * (R̃ \ f)
 end
 function ldiv(Qn::SubQuasiArray{<:Any,2,<:SemiclassicalJacobi,<:Tuple{<:Inclusion,<:Any}}, C::AbstractQuasiArray)
     _,jr = parentindices(Qn)
     Q = parent(Qn)
     if iszero(Q.a) && iszero(Q.b) && Q.c == -one(eltype(Q.t))
-        # TODO: due to a stdlib error this won't work with bigfloat as is
         T = typeof(Q.t)
         R = legendre(zero(T)..one(T))
         B = neg1c_tolegendre(Q.t)
         return (B[jr,jr] \ (R[:,jr] \ C))
     end
-    R = SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1))
-    R̃ = toclassical(R)
+
+    R̃ = toclassical(SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1)))
     (Q \ R̃)[jr,jr] * (R̃[:,jr] \ C)
 end
 

src/neg1c.jl

@@ -4,7 +4,10 @@
 ######
 
 # compute n-th coefficient from direct evaluation formula
-αdirect(n, t::T) where T = 2*n*_₂F₁((one(T)+n)/2,(n+2)/2,(2*n+3)/2,1/t^2)/(t*2*(1+2*n)*_₂F₁(n/2,(n+one(T))/2,(2*n+one(T))/2,1/t^2))
+function αdirect(n, tt::T) where T
+    t = big(tt)
+    return convert(T,2*n*_₂F₁((one(T)+n)/2,(n+2)/2,(2*n+3)/2,1/t^2)/(t*2*(1+2*n)*_₂F₁(n/2,(n+one(T))/2,(2*n+one(T))/2,1/t^2)))
+end
 
 # inital value n=0 for α_{n,n-1}(t) coefficients
 initialα(t) = t-2/(log1p(t)-log(t-1))

src/twoband.jl

@@ -185,15 +185,14 @@ end
 
 function divmul(R::TwoBandJacobi, D::Derivative, P::TwoBandJacobi)
     T = promote_type(eltype(R), eltype(P))
-    ρ=convert(T,R.ρ); t=inv(one(T)-ρ^2)
+    ρ = convert(T,R.ρ); t=inv(one(T)-ρ^2)
     x = axes(R.Q, 1)
 
     Dₑ = -T(2) * t .* (R.Q \ (Derivative(x)*P.P))
     Dₒ = -T(2) .* (HalfWeighted{:c}(R.P) \ (Derivative(x)* HalfWeighted{:c}(P.Q)))
-
-    d = Dₑ.args[1] .* Dₑ.args[2].parent.data
+    d = Dₑ.args[1][1,1] .* Dₑ.args[2].parent.data
     (dₑ, dlₑ) = d[1,:], d[2, :]
-    (dₒ, dlₒ) = Dₒ.data[1,:], Dₒ.data[2,:]
+    (dₒ, dlₒ) = Vcat(-one(T),Dₒ[band(1)]), Dₒ[band(0)]
     BandedMatrix(1=>Interlace(dlₒ, -dₑ), 3=>Interlace(-dₒ[2:∞],dlₑ))
 end
 
@@ -205,7 +204,7 @@ function divmul(R::TwoBandJacobi, D::Derivative, HP::HalfWeighted{:ab,<:Any,<:Tw
     Dₑ = -2*(one(T)-ρ^2) .* (R.Q \ (Derivative(axes(R.Q,1))*HalfWeighted{:ab}(HP.P.P)))
     D₀ = -2*(one(T)-ρ^2)^2 .* (Weighted(R.P) \ (Derivative(axes(R.P,1))*Weighted(HP.P.Q)))
 
-    (dₑ, dlₑ, d₀, dl₀) = Dₑ.data[1,:], Dₑ.data[2,:], D₀.data[1,:], D₀.data[2,:]
+    (dₑ, dlₑ, d₀, dl₀) = Dₑ[band(0)], Dₑ[band(-1)], D₀[band(-1)], D₀[band(-2)]
     BandedMatrix(-1=>Interlace(dₑ, -d₀), -3=>Interlace(-dlₑ, dl₀))
 end
 

test/runtests.jl

@@ -59,6 +59,12 @@ end
     @test w == jacobiexpansion(w)
 end
 
+@testset "Evaluations and Type declaration" begin
+    @test SemiclassicalJacobi(1.1,0,0,4)[0.2, 1] ≈ SemiclassicalJacobi{Float64}(1.1,0,0,4)[0.2, 1]
+    @test SemiclassicalJacobi(1.013,0,0,4)[0.1993, 3] ≈ SemiclassicalJacobi{Float64}(1.013,0,0,4)[0.1993, 3]
+    @test SemiclassicalJacobi(1.013,2,2,2.3)[0.7, 4] ≈ SemiclassicalJacobi{Float64}(1.013,2,2,2.3)[0.7, 4]
+end
+
 @testset "Half-range Chebyshev" begin
     @testset "T and W" begin
         T = SemiclassicalJacobi(2, -1/2, 0, -1/2)
@@ -469,6 +475,39 @@ end
     end
 end
 
+@testset "Weighted Conversion" begin
+    @testset "low c" begin
+        P = SemiclassicalJacobi(1.1,0,0,10)
+        Q = SemiclassicalJacobi(1.1,1,1,10)
+        R = Q \ P
+        L = Weighted(P) \ Weighted(Q)
+        C = R./L'
+        @test C[1,1] ≈ C[5,5]  ≈ C[8,8]
+
+        P = SemiclassicalJacobi(1.6,0,0,10)
+        Q = SemiclassicalJacobi(1.6,2,2,10)
+        R = Q \ P
+        L = Weighted(P) \ Weighted(Q)
+        C = R./L'
+        @test C[1,1] ≈ C[5,5]  ≈ C[8,8]
+
+        P = SemiclassicalJacobi(1.1213,0,0,10)
+        Q = SemiclassicalJacobi(1.1213,2,3,10)
+        R = Q \ P
+        L = Weighted(P) \ Weighted(Q)
+        C = R./L'
+        @test C[1,1] ≈ C[5,5]  ≈ C[8,8]
+    end
+    @testset "high c" begin
+        P = SemiclassicalJacobi(1.1,0,0,100)
+        Q = SemiclassicalJacobi(1.1,1,1,100)
+        R = Q \ P
+        L = Weighted(P) \ Weighted(Q)
+        C = R./L'
+        @test C[1,1] ≈ C[5,5]  ≈ C[8,8]
+    end
+end
+
 @testset "Conversion operators via decomposition" begin
     t = 1.2
     P = SemiclassicalJacobi(t, 1, 0, 1)