Cholesky instead of Lanczos (#131)

* Cholesky instead of Lanczos

* Remove Lanczos

* Get more tests passing

* more tests pass

* cleanup test

* v0.6

* Update runtests.jl

* Update runtests.jl

* Update SemiclassicalOrthogonalPolynomials.jl

* fix whitespace

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.5.11"
+version = "0.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -20,13 +20,13 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 [compat]
 ArrayLayouts = "1"
 BandedMatrices = "0.17, 1"
-ClassicalOrthogonalPolynomials = "0.13, 0.14"
+ClassicalOrthogonalPolynomials = "0.14.3"
 ContinuumArrays = "0.18"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.14, 0.15"
-InfiniteLinearAlgebra = "0.8, 0.9"
-LazyArrays = "2"
+InfiniteLinearAlgebra = "0.9.1"
+LazyArrays = "2.4"
 QuasiArrays = "0.11"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"

src/SemiclassicalOrthogonalPolynomials.jl

@@ -8,19 +8,19 @@ import Base: getindex, axes, size, \, /, *, +, -, summary, show, ==, copy, sum,
 import ArrayLayouts: MemoryLayout, ldiv, diagonaldata, subdiagonaldata, supdiagonaldata
 import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedLayout, _BandedMatrix
 import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, ApplyArray,
-                    AccumulateAbstractVector, LazyVector, AbstractCachedMatrix, BroadcastLayout
+                    AccumulateAbstractVector, LazyVector, AbstractCachedMatrix, BroadcastLayout, simplifiable
 import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedOPLayout, MappedOPLayout,
                                         Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,
                                         OrthogonalPolynomialRatio, Weighted, AbstractWeightLayout, UnionDomain, oneto, WeightedBasis, HalfWeighted,
-                                        golubwelsch, AbstractOPLayout, weight, cholesky_jacobimatrix, qr_jacobimatrix, isnormalized
+                                        golubwelsch, AbstractOPLayout, weight, cholesky_jacobimatrix, qr_jacobimatrix, isnormalized, ConvertedOrthogonalPolynomial, AbstractNormalizedOPLayout
 
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid, equals_layout, ExpansionLayout
 import FillArrays: SquareEye
 import HypergeometricFunctions: _₂F₁general2, _₂F₁
 import SpecialFunctions: beta
 
-export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio
+export Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio
 
 """"
     SemiclassicalJacobiWeight(t, a, b, c)
@@ -131,7 +131,7 @@ WeightedSemiclassicalJacobi{T}(t, a, b, c, P...) where T = SemiclassicalJacobiWe
 
 # Returns α, β such that P1(x) = β(x-α)
 function _linear_coefficients(t, a, b, c)
-    # beta(a + 1, b + 1) * t^c * _₂F₁(a + 1, -c, a + b + 2, 1/t) is the integral ∫₀¹ wᵗ⁽ᵃᵇᶜ⁾(x) dx 
+    # beta(a + 1, b + 1) * t^c * _₂F₁(a + 1, -c, a + b + 2, 1/t) is the integral ∫₀¹ wᵗ⁽ᵃᵇᶜ⁾(x) dx
     Γᵃ = beta(a + 1, b + 1) * _₂F₁(a + 1, -c, a + b + 2, 1/t)
     Γᵃ⁺¹ = beta(a + 2, b + 1) * _₂F₁(a + 2, -c, a + b + 3, 1/t)
     Γᵃ⁺² = beta(a + 3, b + 1) * _₂F₁(a + 3, -c, a + b + 4, 1/t)
@@ -150,7 +150,7 @@ function semiclassical_jacobimatrix(t, a, b, c)
         C = -(N)./(N.*4 .- 2)
         B = Vcat((α[1]^2*3-α[1]*α[2]*2-1)/6 , -(N)./(N.*4 .+ 2).*α[2:end]./α)
         return SymTridiagonal(A, sqrt.(B.*C)) # if J is Tridiagonal(c,a,b) then for norm. OPs it becomes SymTridiagonal(a, sqrt.( b.* c))
-    elseif b == -one(T) 
+    elseif b == -one(T)
         J′ = semiclassical_jacobimatrix(t, a, one(b), c)
         J′a, J′b = diagonaldata(J′), supdiagonaldata(J′)
         A = Vcat(one(T), J′a[1:end])
@@ -168,7 +168,7 @@ function semiclassical_jacobimatrix(t, a, b, c)
             return SemiclassicalJacobi.(t, a, b, 0:Int(c))[end].X
         else # if c is not an integer, use Lanczos
             x = axes(P,1)
-            return jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), jacobi(b, a, UnitInterval{T}())))
+            return cholesky_jacobimatrix(@.(x^a * (1-x)^b * (t-x)^c), P)
         end
     end
 end
@@ -178,11 +178,11 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
     Δb = b-Q.b
     Δc = c-Q.c
 
-    # special cases 
+    # special cases
     if iszero(a) && iszero(b) && c == -one(eltype(Q.t)) # (a,b,c) = (0,0,-1) special case
         return semiclassical_jacobimatrix(Q.t, zero(Q.t), zero(Q.t), c)
     elseif iszero(Δa) && iszero(Δc) && Δb == 2 && b == 1
-        # When going from P[t, a, -1, c] to P[t, a, 1, c], you can just take 
+        # When going from P[t, a, -1, c] to P[t, a, 1, c], you can just take
         return SymTridiagonal(Q.X.d[2:end], Q.X.du[2:end])
     elseif iszero(c) # classical Jacobi polynomial special case
         return jacobimatrix(Normalized(jacobi(b, a, UnitInterval{eltype(Q.t)}())))
@@ -233,15 +233,14 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
     end
 end
 
-LanczosPolynomial(P::SemiclassicalJacobi{T}) where T =
-    LanczosPolynomial(orthogonalityweight(P), Normalized(jacobi(P.b, P.a, UnitInterval{T}())), P.X.dv.data)
+ConvertedOrthogonalPolynomial(P::SemiclassicalJacobi{T}) where T = ConvertedOrthogonalPolynomial(orthogonalityweight(P), P.X, parent(P.X.dv).U, parent(P.X.dv).P)
 
 """
     toclassical(P::SemiclassicalJacobi)
 
-gives either a mapped `Jacobi` or `LanczosPolynomial` version of `P`.
+gives either a mapped `Jacobi` or `CholeskyPolynomial` version of `P`.
 """
-toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : LanczosPolynomial(P)
+toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : ConvertedOrthogonalPolynomial(P)
 
 copy(P::SemiclassicalJacobi) = P
 axes(P::SemiclassicalJacobi{T}) where T = (Inclusion(UnitInterval{T}()),OneToInf())
@@ -281,7 +280,7 @@ function op_lowering(Q, y)
 end
 
 function semijacobi_ldiv(Q, P::SemiclassicalJacobi)
-    if P.a ≤ 0 && P.b ≤ 0 && P.c ≤ 0
+    if P.a ≤ 0 && P.b ≤ 0 && P.c ≤ 0
         P̃ = toclassical(P)
         (Q \ P̃)/_p0(P̃)
     else
@@ -298,7 +297,7 @@ end
 """
     semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
 
-Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q` in a single step via decomposition. 
+Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q` in a single step via decomposition.
 Numerically unstable if the parameter modification is large. Typically one should instead use `P \\ Q` which is equivalent to `semijacobi_ldiv(P,Q)` and proceeds step by step.
 """
 function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
@@ -365,7 +364,7 @@ end
 """
     semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
 
-Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q`. Integer distances are covered by decomposition methods, for non-integer cases a Lanczos fallback is attempted.
+Returns conversion operator from SemiclassicalJacobi `P` to SemiclassicalJacobi `Q`. Integer distances are covered by decomposition methods.
 """
 function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
     @assert Q.t ≈ P.t
@@ -396,7 +395,7 @@ function semijacobi_ldiv(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
             QQ = SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1+iseven(Δc), P)
             return ApplyArray(*,semijacobi_ldiv_direct(Q, QQ),semijacobi_ldiv(QQ, P))
         end
-    else # fallback to Lancos
+    else # fallback
         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)
@@ -406,8 +405,8 @@ end
 struct SemiclassicalJacobiLayout <: AbstractOPLayout end
 MemoryLayout(::Type{<:SemiclassicalJacobi}) = SemiclassicalJacobiLayout()
 
-copy(L::Ldiv{<:NormalizedOPLayout,SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}(L.A, L.B))
-copy(L::Ldiv{SemiclassicalJacobiLayout,<:NormalizedOPLayout}) = copy(Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
+copy(L::Ldiv{<:AbstractNormalizedOPLayout,SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}(L.A, L.B))
+copy(L::Ldiv{SemiclassicalJacobiLayout,<:AbstractNormalizedOPLayout}) = copy(Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
 
 copy(L::Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},BasisLayout}(L.A, L.B))
 copy(L::Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
@@ -419,6 +418,7 @@ copy(L::Ldiv{WeightedOPLayout,SemiclassicalJacobiLayout}) = copy(Ldiv{WeightedOP
 copy(L::Ldiv{SemiclassicalJacobiLayout,WeightedOPLayout}) = copy(Ldiv{BasisLayout,WeightedOPLayout}(L.A, L.B))
 
 
+simplifiable(L::Ldiv{SemiclassicalJacobiLayout,<:AbstractBasisLayout}) = Val(true)
 copy(L::Ldiv{SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
 copy(L::Ldiv{SemiclassicalJacobiLayout,<:AbstractBasisLayout}) = semijacobi_ldiv(L.A, L.B)
 copy(L::Ldiv{SemiclassicalJacobiLayout,BroadcastLayout{typeof(*)}}) = semijacobi_ldiv(L.A, L.B)
@@ -435,25 +435,25 @@ function copy(L::Ldiv{SemiclassicalJacobiLayout,SemiclassicalJacobiLayout})
     (inv(M_Q) * L') * M_P
 end
 
-function \(A::SemiclassicalJacobi, B::SemiclassicalJacobi{T}) where {T} 
+function \(A::SemiclassicalJacobi, B::SemiclassicalJacobi{T}) where {T}
     if A.b == -1 && B.b ≠ -1
-        return UpperTriangular(ApplyArray(inv, B \ A)) 
+        return UpperTriangular(ApplyArray(inv, B \ A))
     elseif B.b == -1 && A.b ≠ -1
         # First convert Bᵗᵃ⁻¹ᶜ into Bᵗᵃ⁰ᶜ
-        Bᵗᵃ⁰ᶜ = SemiclassicalJacobi(B.t, B.a, zero(B.b), B.c, A) 
+        Bᵗᵃ⁰ᶜ = SemiclassicalJacobi(B.t, B.a, zero(B.b), B.c, A)
         Bᵗᵃ¹ᶜ = SemiclassicalJacobi(B.t, B.a, one(B.a), B.c, A)
         Rᵦₐ₁ᵪᵗᵃ⁰ᶜ = Weighted(Bᵗᵃ⁰ᶜ) \ Weighted(Bᵗᵃ¹ᶜ)
         b1 = Rᵦₐ₁ᵪᵗᵃ⁰ᶜ[band(0)]
         b0 = Vcat(one(T), Rᵦₐ₁ᵪᵗᵃ⁰ᶜ[band(-1)])
         Rᵦₐ₋₁ᵪᵗᵃ⁰ᶜ = Bidiagonal(b0, b1, :U)
-        # Then convert Bᵗᵃ⁰ᶜ into A and complete 
+        # Then convert Bᵗᵃ⁰ᶜ into A and complete
         Rₐ₀ᵪᴬ = UpperTriangular(A \ Bᵗᵃ⁰ᶜ)
         return ApplyArray(*, Rₐ₀ᵪᴬ, Rᵦₐ₋₁ᵪᵗᵃ⁰ᶜ)
     elseif A.b == B.b == -1
         Bᵗᵃ¹ᶜ = SemiclassicalJacobi(B.t, B.a, one(B.b), B.c, B)
         Aᵗᵃ¹ᶜ = SemiclassicalJacobi(A.t, A.a, one(A.b), A.c, A)
         Rₐ₁ᵪᵗᵘ¹ᵛ = Aᵗᵃ¹ᶜ \ Bᵗᵃ¹ᶜ
-        # Make 1 ⊕ Rₐ₁ᵪᵗᵘ¹ᵛ 
+        # Make 1 ⊕ Rₐ₁ᵪᵗᵘ¹ᵛ
         V = eltype(Rₐ₁ᵪᵗᵘ¹ᵛ)
         Rₐ₋₁ᵪᵗᵘ⁻¹ᵛ = Vcat(
             Hcat(one(V), Zeros{V}(1, ∞)),
@@ -464,8 +464,8 @@ function \(A::SemiclassicalJacobi, B::SemiclassicalJacobi{T}) where {T}
         return semijacobi_ldiv(A, B)
     end
 end
-\(A::LanczosPolynomial, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
-\(A::SemiclassicalJacobi, B::LanczosPolynomial) = semijacobi_ldiv(A, B)
+\(A::ConvertedOrthogonalPolynomial, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
+\(A::SemiclassicalJacobi, B::ConvertedOrthogonalPolynomial) = semijacobi_ldiv(A, B)
 function \(w_A::WeightedSemiclassicalJacobi{T}, w_B::WeightedSemiclassicalJacobi{T}) where T
     wA,A = w_A.args
     wB,B = w_B.args
@@ -479,12 +479,12 @@ function \(w_A::WeightedSemiclassicalJacobi{T}, w_B::WeightedSemiclassicalJacobi
         @assert A.a + 1 == B.a && A.c + 1 == B.c
         Q = SemiclassicalJacobi(B.t, B.a, one(B.b), B.c, B)
         P = SemiclassicalJacobi(A.t, A.a, one(A.b), A.c, A)
-        wP = Weighted(P) 
+        wP = Weighted(P)
         wQ = Weighted(Q)
         R22 = wP \ wQ
         α, β = _linear_coefficients(P.t, P.a, P.b, P.c)
-        ℓ₁ = A.t - 1 
-        ℓ₂ = 1 + α - A.t 
+        ℓ₁ = A.t - 1
+        ℓ₂ = 1 + α - A.t
         ℓ₃ = inv(β)
         d0 = Vcat(ℓ₁, R22[band(0)])
         d1 = Vcat(ℓ₂, R22[band(-1)])
@@ -505,13 +505,13 @@ function \(w_A::WeightedSemiclassicalJacobi{T}, w_B::WeightedSemiclassicalJacobi
     elseif wA.a == wB.a && wA.b == wB.b && wA.c+1 == wB.c
         @assert A.a == B.a && A.b == B.b && A.c+1 == B.c
         -op_lowering(A,A.t)
-    elseif wA.a+1 ≤ wB.a
+    elseif wA.a+1 ≤ wB.a
         C = SemiclassicalJacobi(A.t, A.a+1, A.b, A.c, A)
         w_C = SemiclassicalJacobiWeight(wA.t, wA.a+1, wA.b, wA.c) .* C
         L_2 = w_C \ w_B
         L_1 = w_A \ w_C
         L_1 * L_2
-    elseif wA.b+1 ≤ wB.b
+    elseif wA.b+1 ≤ wB.b
         C = SemiclassicalJacobi(A.t, A.a, A.b+1, A.c, A)
         w_C = SemiclassicalJacobiWeight(wA.t, wA.a, wA.b+1, wA.c) .* C
         L_2 = w_C \ w_B
@@ -529,19 +529,19 @@ end
 \(w_A::Weighted{<:Any,<:SemiclassicalJacobi}, w_B::Weighted{<:Any,<:SemiclassicalJacobi}) = convert(WeightedBasis, w_A) \ convert(WeightedBasis, w_B)
 
 function \(w_A::HalfWeighted{lr, T, <:SemiclassicalJacobi}, B::AbstractQuasiArray{V}) where {lr, T, V}
-    WP = convert(WeightedBasis, w_A) 
-    w_A.P.b  ≠ -1 && return WP \ B # no need to special case here 
+    WP = convert(WeightedBasis, w_A)
+    w_A.P.b  ≠ -1 && return WP \ B # no need to special case here
     !iszero(WP.args[1].b) && throw(ArgumentError("Cannot expand in a weighted basis including 1/(1-x)."))
-    # To expand f(x) = w(x)P(x)𝐟, note that P = [1 (1-x)Q] so 
+    # To expand f(x) = w(x)P(x)𝐟, note that P = [1 (1-x)Q] so
     #   f(x) = w(x)[1 (1-x)Q(x)][f₀; 𝐟₁] = w(x)f₀ + w(x)(1-x)Q(x)𝐟₁. Thus,
-    #   f(1) = w(1)f₀ ⟹ f₀ = f(1) / w(1) 
-    #   Then, f(x) - w(x)f₀ = w(x)(1-x)Q(x)𝐟₁, so that 𝐟₁ is just the expansion of 
+    #   f(1) = w(1)f₀ ⟹ f₀ = f(1) / w(1)
+    #   Then, f(x) - w(x)f₀ = w(x)(1-x)Q(x)𝐟₁, so that 𝐟₁ is just the expansion of
     #   f(x) - w(x)f₀ in the w(x)(1-x)Q(x) basis.
-    w, P = WP.args 
-    f₀ = B[end] / w[end] 
+    w, P = WP.args
+    f₀ = B[end] / w[end]
     C = B - w * f₀
     Q = SemiclassicalJacobiWeight(w.t, w.a, one(w.b), w.c) .* SemiclassicalJacobi(P.t, P.a, one(P.b), P.c, P)
-    f = Q \ C 
+    f = Q \ C
     return Vcat(f₀, f)
 end
 
@@ -563,8 +563,8 @@ function ldiv(Q::SemiclassicalJacobi, f::AbstractQuasiVector)
     elseif isinteger(Q.a) && (isinteger(Q.b) && Q.b ≥ 0) && 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)
-    elseif isinteger(Q.a) && isone(-Q.b) && isinteger(Q.c) 
-        return semijacobi_ldiv(Q, f) # jacobi(< 0, Q.a) fails in the method above. jacobi(-1, 0) also leads to NaNs in coefficients 
+    elseif isinteger(Q.a) && isone(-Q.b) && isinteger(Q.c)
+        return semijacobi_ldiv(Q, f) # jacobi(< 0, Q.a) fails in the method above. jacobi(-1, 0) also leads to NaNs in coefficients
     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)
@@ -639,7 +639,7 @@ function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
 end
 
 function _getsecondifpossible(v)
-    length(v) > 1 && return v[2] 
+    length(v) > 1 && return v[2]
     return v[1]
 end
 
@@ -651,9 +651,9 @@ end
 
 Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
 Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
-Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) = 
+Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) =
     SemiclassicalJacobiFamily(t, a, b, c)
-Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) where T = 
+Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Union{AbstractUnitRange,Number}, b::Union{AbstractUnitRange,Number}, c::Union{AbstractUnitRange,Number}) where T =
     SemiclassicalJacobiFamily{T}(t, a, b, c)
 
 
@@ -664,7 +664,7 @@ function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::Abstract
     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 
+        # 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])
         P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), Pprev)
@@ -704,7 +704,7 @@ function SemiclassicalJacobiCWeightFamily{T}(t::Number, a::Number, b::Number, c:
     return SemiclassicalJacobiCWeightFamily{T}(SemiclassicalJacobiWeight.(t,a:a,b:b,c), t, a, b, c)
 end
 
-Base.broadcasted(::Type{SemiclassicalJacobiWeight}, t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) = 
+Base.broadcasted(::Type{SemiclassicalJacobiWeight}, t::Number, a::Number, b::Number, c::Union{AbstractUnitRange,Number}) =
 SemiclassicalJacobiCWeightFamily(t, a, b, c)
 
 _unweightedsemiclassicalsum(a,b,c,t) = pFq((a+1,-c),(a+b+2, ), 1/t)
@@ -714,7 +714,7 @@ function Base.broadcasted(::typeof(sum), W::SemiclassicalJacobiCWeightFamily{T})
     cmin = minimum(c); cmax = maximum(c);
     @assert isinteger(cmax) && isinteger(cmin)
     # This is needed at high parameter values.
-    # Manually setting setprecision(2048) allows accurate computation even for very high c. 
+    # Manually setting setprecision(2048) allows accurate computation even for very high c.
     t,a,b = convert(BigFloat,t),convert(BigFloat,a),convert(BigFloat,b)
     F = zeros(BigFloat,cmax+1)
     F[1] = _unweightedsemiclassicalsum(a,b,0,t) # c=0

test/runtests.jl

@@ -72,6 +72,10 @@ end
     R = Q \ P
     Ralt = semijacobi_ldiv_direct(Q,P)
     @test R[1:20,1:20] ≈ Ralt[1:20,1:20]
+
+    @test (P \ legendre(0..1))[1:10,1:10] * (legendre(0..1) \ Normalized(legendre(0..1)))[1:10,1:10] ≈ (P \ Normalized(legendre(0..1)))[1:10,1:10]
+    @test_broken (Normalized(legendre(0..1)) \ legendre(0..1))[1:10,1:10] * (legendre(0..1) \ P)[1:10,1:10] ≈ (Normalized(legendre(0..1)) \ P)[1:10,1:10]
+
     # set 2
     P = SemiclassicalJacobi(1.23,4,1,2)
     Q = SemiclassicalJacobi(1.23,7,4,6)
@@ -274,7 +278,7 @@ end
         U = SemiclassicalJacobi(2, 1/2, 0, 1/2, T)
         x = axes(T,1)
 
-        ucfs = (T \ exp.(x))
+        ucfs = T \ exp.(x)
         u = T * ucfs
         @test u[0.1] ≈ exp(0.1)
         @test T[:,1:20] \ exp.(x) ≈ ucfs[1:20]