Computing connection matrix for decrementing returns Nothing

[DanielVandH]

julia> using SemiclassicalOrthogonalPolynomials

julia> t, a, b, c = 2.0, 1.0, -1.0, 1.0
(2.0, 1.0, -1.0, 1.0)

julia> A = SemiclassicalJacobi(t, a + 1, 0, c + 1)
SemiclassicalJacobi{Float64}(2.0, 2.0, 0.0, 2.0, [0.6874999999999998 0.21260501606769022 …  … 0.21260501606769022 0.5717592592592593 …  … … ; … ; ])

julia> B = SemiclassicalJacobi(t, a + 1, 1, c + 1)
SemiclassicalJacobi{Float64}(2.0, 2.0, 1.0, 2.0, [0.542857142857143 0.2025349554108261 …  … 0.2025349554108261 0.5284529732290925 …  … … ; … ; ])

julia> A \ B

The issue is with semijacobi_ldiv:

function semijacobi_ldiv(Qt, Qa, Qb, Qc, P::SemiclassicalJacobi)
    @assert Qt ≈ P.t
    (Qt ≈ P.t) && (Qa ≈ P.a) && (Qb ≈ P.b) && (Qc ≈ P.c) && return SymTridiagonal(Ones(∞),Zeros(∞))
    Δa = Qa-P.a
    Δb = Qb-P.b
    Δc = Qc-P.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(Qt, Qa, Qb, Qc, P)
        elseif Δa > 0  # iterative modification by 1
            M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa-1-iseven(Δa), Qb, Qc, P)),semijacobi_ldiv(Qt, Qa-1-iseven(Δa), Qb, Qc, P))
        elseif Δb > 0
            M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb-1-iseven(Δb), Qc, P)),semijacobi_ldiv(Qt, Qa, Qb-1-iseven(Δb), Qc, P))
        elseif Δc > 0
            M = ApplyArray(*,semijacobi_ldiv_direct(Qt, Qa, Qb, Qc, SemiclassicalJacobi(Qt, Qa, Qb, Qc-1-iseven(Δc), P)),semijacobi_ldiv(Qt, Qa, Qb, Qc-1-iseven(Δc), P))
        end
    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̃ \ SemiclassicalJacobi(Qt, Qa, Qb, Qc))
    end
end

Since Δb < 0, none of the conditions in the first block of if statements are satisfied, but since it ends in an elseif it just returns Nothing. I'll have to try and make it handle negatives one at a time and then inverse at the end, just like is done with e.g. Jacobi.

[TSGut]

Lowering used to be implemented before we moved the package to decomposition based connection coefficients. You will see a TODO comment in semiclassical_jacobimatrix about this as well. It's in principle not an issue to implement this using QL or reverse Cholesky, would be nice to have.

[DanielVandH]

Missed that TODO comment. I think for now I'm going to just look at implementing the inv approach into ldiv rather than working with QL / reverse Cholesky itself, mimicing what's done in Jacobi

[DanielVandH]

This works now on master