diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 10f85a5..76ebcc3 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -174,6 +174,11 @@ function semiclassical_jacobimatrix(t, a, b, c) end end +function semiclassical_jacobimatrix_raise_c_by_2(Q) + X = jacobimatrix(Q) + qr_jacobimatrix(Q.t*I-X,X)[1] +end + function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c) Δa = a-Q.a Δb = b-Q.b @@ -202,7 +207,7 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c) elseif iszero(Δa) && isone(Δb/2) && iszero(Δc) qr_jacobimatrix(I-X,X)[1] elseif iszero(Δa) && iszero(Δb) && isone(Δc/2) - qr_jacobimatrix(Q.t*I-X,X)[1] + semiclassical_jacobimatrix_raise_c_by_2(Q) elseif isone(Δa) && iszero(Δb) && iszero(Δc) # raising by 1 cholesky_jacobimatrix(X,X)[1] elseif iszero(Δa) && isone(Δb) && iszero(Δc) diff --git a/src/family.jl b/src/family.jl index 49c724d..8cdd288 100644 --- a/src/family.jl +++ b/src/family.jl @@ -53,13 +53,21 @@ Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUni _broadcast_getindex(a,k) = a[k] _broadcast_getindex(a::Number,k) = a -function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange) +function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily{T,<:Number,<:Number,<:AbstractUnitRange}, inds::AbstractUnitRange) where T + t,a,b,c = P.t,P.a,P.b,P.c + for k in inds + Pprev = P.data[k-2] + P.data[k] = SemiclassicalJacobi{T}(Pprev.t, Pprev.a, Pprev.b, Pprev.c+2, semiclassical_jacobimatrix_raise_c_by_2(Pprev)) + end + P +end + +function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily{<:Number,<:Number,<:AbstractUnitRange}, inds::AbstractUnitRange) t,a,b,c = P.t,P.a,P.b,P.c - isrange = P.b isa AbstractUnitRange for k in inds # If P.data[k-2] is not normalised (aka b = -1), cholesky fails. With the current design, this is only a problem if P.b # is a range since we can translate between polynomials that both have b = -1. - Pprev = (isrange && P.b[k-2] == -1) ? P.data[k-1] : P.data[k-2] # isrange && P.b[k-2] == -1 could also be !isnormalized(P.data[k-2]) + Pprev = P.b[k-2] == -1 ? P.data[k-1] : P.data[k-2] # isrange && P.b[k-2] == -1 could also be !isnormalized(P.data[k-2]) P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), Pprev) end P