Improvements to Normalized (#8)

* Update normalized.jl

* custom normalized

* Actually super simple... duh...

* Simplify recurrencecoefficients

* Switch to Infinities.jl

* Update Project.toml

* Move to new Bi/Tri/SymTridiagonal

* fix some tests

* tests pass

* Update Project.toml

* Update Project.toml

* Update Project.toml

* Update Project.toml

* Update Project.toml

* Realised Hermite and Stieltjes tests weren't being run... need to come back later and fix bugs

* tests pass

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,35 +18,37 @@ InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "0.5.3"
-BandedMatrices = "0.16.4"
+ArrayLayouts = "0.6.2"
+BandedMatrices = "0.16.5"
 BlockArrays = "0.14.1"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.5"
+ContinuumArrays = "0.6"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"
-InfiniteArrays = "0.9"
-InfiniteLinearAlgebra = "0.4.6"
-IntervalSets = "0.3.1, 0.4, 0.5"
-LazyArrays = "0.20.5"
-QuasiArrays = "0.4.5"
+InfiniteArrays = "0.10"
+InfiniteLinearAlgebra = "0.5.2"
+IntervalSets = "0.5"
+LazyArrays = "0.20.6"
+LazyBandedMatrices = "0.5"
+QuasiArrays = "0.4.6"
 SpecialFunctions = "0.10, 1"
 julia = "1.5"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "LazyBandedMatrices"]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices"]

src/ClassicalOrthogonalPolynomials.jl

@@ -1,7 +1,8 @@
 module ClassicalOrthogonalPolynomials
 using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
     IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
-    InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW
+    InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
+    LazyBandedMatrices
 
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
@@ -10,8 +11,11 @@ import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
-                _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout, AbstractCachedVector
-import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata
+                _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout, 
+                AbstractCachedVector, AbstractCachedMatrix
+import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
+                        subdiagonaldata, diagonaldata, supdiagonaldata
+import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal
 import LinearAlgebra: pinv, factorize, qr, adjoint, transpose
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
 import FillArrays: AbstractFill, getindex_value
@@ -22,7 +26,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, _factorize
 
-import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange
+import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
     AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
@@ -53,7 +57,7 @@ include("interlace.jl")
 cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
 cardinality(::EuclideanDomain) = ℵ₁
 
-transform_ldiv(A, f, ::Tuple{<:Any,Infinity})  = adaptivetransform_ldiv(A, f)
+transform_ldiv(A, f, ::Tuple{<:Any,InfiniteCardinal{0}})  = adaptivetransform_ldiv(A, f)
 
 function chop!(c::AbstractVector, tol::Real)
     @assert tol >= 0
@@ -213,6 +217,17 @@ function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::Sub
     P[kr, :] * view(X,:,jr)
 end
 
+function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
+    T = promote_type(eltype(x), eltype(C))
+    x == axes(C,1) || throw(DimensionMismatch())
+    P = parent(C)
+    kr,_ = parentindices(C)
+    y = axes(P,1)
+    Y = P \ (y .* P)
+    X = kr.A \ (Y     - kr.b * Eye{T}(∞))
+    P[kr, :] * X
+end
+
 function jacobimatrix(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) where T
     P = parent(C)
     kr,jr = parentindices(C)
@@ -250,6 +265,11 @@ function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:
     (parent(A) \ parent(B))[jA, jB]
 end
 
+function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}})
+    axes(A,1) == axes(B,1) || throw(DimensionMismatch())
+    parent(A) \ parent(B)
+end
+
 function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
     w_A,A = arguments(wA)
     w_B,B = arguments(wB)

src/chebyshev.jl

@@ -106,15 +106,11 @@ factorize(L::SubQuasiArray{T,2,<:ChebyshevU,<:Tuple{<:Inclusion,<:OneTo}}) where
 # Jacobi Matrix
 ########
 
-jacobimatrix(C::ChebyshevT{T}) where T =
-    _BandedMatrix(Vcat(Fill(one(T)/2,1,∞),
-                        Zeros{T}(1,∞),
-                        Hcat(one(T), Fill(one(T)/2,1,∞))), ∞, 1, 1)
+jacobimatrix(C::ChebyshevT{T}) where T = 
+    Tridiagonal(Vcat(one(T), Fill(one(T)/2,∞)), Zeros{T}(∞), Fill(one(T)/2,∞))
 
 jacobimatrix(C::ChebyshevU{T}) where T =
-    _BandedMatrix(Vcat(Fill(one(T)/2,1,∞),
-                        Zeros{T}(1,∞),
-                        Fill(one(T)/2,1,∞)), ∞, 1, 1)
+    SymTridiagonal(Zeros{T}(∞), Fill(one(T)/2,∞))
 
 
 
@@ -151,7 +147,7 @@ end
 # Ultraspherical(1)\(D*Chebyshev())
 @simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::ChebyshevT)
     T = promote_type(eltype(D),eltype(S))
-    A = _BandedMatrix((zero(T):∞)', ∞, -1,1)
+    A = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
     ApplyQuasiMatrix(*, ChebyshevU{T}(), A)
 end
 
@@ -163,14 +159,14 @@ function \(U::ChebyshevU, C::ChebyshevT)
     T = promote_type(eltype(U), eltype(C))
     _BandedMatrix(Vcat(-Ones{T}(1,∞)/2,
                         Zeros{T}(1,∞),
-                        Hcat(Ones{T}(1,1),Ones{T}(1,∞)/2)), ∞, 0,2)
+                        Hcat(Ones{T}(1,1),Ones{T}(1,∞)/2)), ℵ₀, 0,2)
 end
 
 function \(w_A::WeightedChebyshevT, w_B::WeightedChebyshevU)
     wA,A = w_A.args
     wB,B = w_B.args
     T = promote_type(eltype(w_A), eltype(w_B))
-    _BandedMatrix(Vcat(Fill(one(T)/2, 1, ∞), Zeros{T}(1, ∞), Fill(-one(T)/2, 1, ∞)), ∞, 2, 0)
+    _BandedMatrix(Vcat(Fill(one(T)/2, 1, ∞), Zeros{T}(1, ∞), Fill(-one(T)/2, 1, ∞)), ℵ₀, 2, 0)
 end
 
 \(T::ChebyshevT, U::ChebyshevU) = inv(U \ T)

src/fourier.jl

@@ -64,7 +64,8 @@ end
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVector{<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
     axes(c,1) == axes(F,1) || throw(DimensionMismatch())
     T = promote_type(eltype(c), eltype(F))
-    F*mortar(Tridiagonal(Vcat([reshape([0; one(T)],2,1)], Fill(Matrix(one(T)/2*I,2,2),∞)),
+    # Use LinearAlgebra.Tridiagonal for now since broadcasting support not complete for LazyBandedMatrices.Tridiagonal
+    F*mortar(LinearAlgebra.Tridiagonal(Vcat([reshape([0; one(T)],2,1)], Fill(Matrix(one(T)/2*I,2,2),∞)),
                         Vcat([zeros(T,1,1)], Fill(Matrix(zero(T)I,2,2),∞)),
                         Vcat([[0 one(T)/2]], Fill(Matrix(one(T)/2*I,2,2),∞))))
 end
@@ -73,7 +74,8 @@ end
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), s::BroadcastQuasiVector{<:Any,typeof(sin),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
     axes(s,1) == axes(F,1) || throw(DimensionMismatch())
     T = promote_type(eltype(s), eltype(F))
-    F*mortar(Tridiagonal(Vcat([reshape([one(T); 0],2,1)], Fill([0 one(T)/2; -one(T)/2 0],∞)),
+    # Use LinearAlgebra.Tridiagonal for now since broadcasting support not complete for LazyBandedMatrices.Tridiagonal
+    F*mortar(LinearAlgebra.Tridiagonal(Vcat([reshape([one(T); 0],2,1)], Fill([0 one(T)/2; -one(T)/2 0],∞)),
                         Vcat([zeros(T,1,1)], Fill(Matrix(zero(T)*I,2,2),∞)),
                         Vcat([[one(T)/2 0]], Fill([0 -one(T)/2; one(T)/2 0],∞))))
 end

src/hermite.jl

@@ -17,10 +17,8 @@ axes(::Hermite{T}) where T = (Inclusion(ℝ), oneto(∞))
 # H_{n+1} = 2x H_n - 2n H_{n-1}
 # 1/2 * H_{n+1} + n H_{n-1} = x H_n 
 # x*[H_0 H_1 H_2 …] = [H_0 H_1 H_2 …] * [0    1; 1/2  0     2; 1/2   0  3; …]   
-function jacobimatrix(H::Hermite{T}) where T
-    # X = BandedMatrix(1 => 1:∞, -1 => Fill(one(T)/2,∞))
-    _BandedMatrix(Vcat((0:∞)', Zeros(1,∞), Fill(one(T)/2,1,∞)), ∞, 1, 1)
-end
+jacobimatrix(H::Hermite{T}) where T = Tridiagonal(Fill(one(T)/2,∞), Zeros{T}(∞), one(T):∞)
+recurrencecoefficients(H::Hermite{T}) where T = Fill{T}(2,∞), Zeros{T}(∞), zero(T):2:∞
 
 @simplify function *(Ac::QuasiAdjoint{<:Any,<:Hermite}, B::WeightedBasis{<:Any,<:HermiteWeight,<:Hermite})  
     T = promote_type(eltype(Ac), eltype(B))
@@ -33,6 +31,6 @@ end
 
 @simplify function *(D::Derivative, H::Hermite)
     T = promote_type(eltype(D),eltype(H))
-    D = _BandedMatrix((zero(T):2:∞)', ∞, -1,1)
+    D = _BandedMatrix((zero(T):2:∞)', ℵ₀, -1,1)
     H*D
 end

src/jacobi.jl

@@ -213,7 +213,7 @@ end
 # Jacobi Matrix
 ########
 
-jacobimatrix(::Legendre{T}) where T = _BandedMatrix(Vcat(((zero(T):∞)./(1:2:∞))', Zeros{T}(1,∞), ((one(T):∞)./(1:2:∞))'), ∞, 1,1)
+jacobimatrix(::Legendre{T}) where T =  Tridiagonal((one(T):∞)./(1:2:∞), Zeros{T}(∞), (one(T):∞)./(3:2:∞))
 
 # These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
 function recurrencecoefficients(::Legendre{T}) where T
@@ -226,9 +226,9 @@ function jacobimatrix(J::Jacobi)
     n = 0:∞
     B = Vcat(2 / (a+b+2),  2 .* (n .+ 2) .* (n .+ (a+b+2)) ./ ((2n .+ (a+b+3)) .* (2n .+ (a+b+4))))
     A = Vcat((b-a) / (a+b+2), (b^2-a^2) ./ ((2n .+ (a+b+2)) .* (2n .+ (a+b+4))))
-    C = 2 .* (n .+ a) .* (n .+ b) ./ ((2n .+ (a+b)) .* (2n .+ (a+b+1)))
+    C = 2 .* (n .+ (a + 1)) .* (n .+ (b + 1)) ./ ((2n .+ (a+b+2)) .* (2n .+ (a+b+3)))
 
-    _BandedMatrix(Vcat(C',A',B'), ∞, 1,1)
+    Tridiagonal(B,A,C)
 end
 
 function recurrencecoefficients(P::Jacobi)
@@ -260,17 +260,13 @@ function \(A::Jacobi, B::Jacobi)
     if A.a ≈ a && A.b ≈ b
         Eye{T}(∞)
     elseif isone(-a-b) && A.a == a && A.b == b+1
-        _BandedMatrix(Vcat((((0:∞) .+ a)./((1:2:∞) .+ (a+b)))',
-                            Vcat(1,((2:∞) .+ (a+b))./((3:2:∞) .+ (a+b)))'), ∞, 0,1)
+        Bidiagonal(Vcat(1, ((2:∞) .+ (a+b)) ./ ((3:2:∞) .+ (a+b))), ((1:∞) .+ a) ./ ((3:2:∞) .+ (a+b)), :U)
     elseif isone(-a-b) && A.a == a+1 && A.b == b
-        _BandedMatrix(Vcat((-((0:∞) .+ b)./((1:2:∞) .+ (a+b)))',
-                            Vcat(1,((2:∞) .+ (a+b))./((3:2:∞) .+ (a+b)))'), ∞, 0,1)
+        Bidiagonal(Vcat(1, ((2:∞) .+ (a+b)) ./ ((3:2:∞) .+ (a+b))), -((1:∞) .+ b) ./ ((3:2:∞) .+ (a+b)), :U)
     elseif A.a == a && A.b == b+1
-        _BandedMatrix(Vcat((((0:∞) .+ a)./((1:2:∞) .+ (a+b)))',
-                            (((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)))'), ∞, 0,1)
+        Bidiagonal(((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)), ((1:∞) .+ a)./((3:2:∞) .+ (a+b)), :U)
     elseif A.a == a+1 && A.b == b
-        _BandedMatrix(Vcat((-((0:∞) .+ b)./((1:2:∞) .+ (a+b)))',
-                            (((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)))'), ∞, 0,1)
+        Bidiagonal(((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)), -((1:∞) .+ b)./((3:2:∞) .+ (a+b)), :U)
     elseif A.a ≥ a+1
         J = Jacobi(a+1,b)
         (A \ J) * (J \ B)
@@ -310,9 +306,9 @@ function \(w_A::WeightedJacobi, w_B::WeightedJacobi)
     if wA == wB
         A \ B
     elseif B.a == A.a && B.b == A.b+1 && wB.b == wA.b+1 && wB.a == wA.a
-        _BandedMatrix(Vcat((((2:2:∞) .+ 2A.b)./((2:2:∞) .+ (A.a+A.b)))', ((2:2:∞)./((2:2:∞) .+ (A.a+A.b)))'), ∞, 1,0)
+        Bidiagonal(((2:2:∞) .+ 2A.b)./((2:2:∞) .+ (A.a+A.b)), (2:2:∞)./((2:2:∞) .+ (A.a+A.b)), :L)
     elseif B.a == A.a+1 && B.b == A.b && wB.b == wA.b && wB.a == wA.a+1
-        _BandedMatrix(Vcat((((2:2:∞) .+ 2A.a)./((2:2:∞) .+ (A.a+A.b)))', -((2:2:∞)./((2:2:∞) .+ (A.a+A.b)))'), ∞, 1,0)
+        Bidiagonal(((2:2:∞) .+ 2A.a)./((2:2:∞) .+ (A.a+A.b)), -(2:2:∞)./((2:2:∞) .+ (A.a+A.b)), :L)
     elseif wB.a ≥ wA.a+1
         J = JacobiWeight(wB.a-1,wB.b) .* Jacobi(B.a-1,B.b)
         (w_A\J) * (J\w_B)
@@ -332,7 +328,7 @@ end
 
 # Jacobi(a+1,b+1)\(D*Jacobi(a,b))
 @simplify function *(D::Derivative{<:Any,<:AbstractInterval}, S::Jacobi)
-    A = _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ∞, -1,1)
+    A = _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ℵ₀, -1,1)
     ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), A)
 end
 
@@ -343,13 +339,13 @@ end
     if w.a == 0 && w.b == 0
         D*S
     elseif iszero(w.a) && w.b == b #L_6
-        A = _BandedMatrix((b:∞)', ∞, 0,0)
+        A = Diagonal(b:∞)
         ApplyQuasiMatrix(*, JacobiWeight(w.a,b-1) .* Jacobi(a+1,b-1), A)
     elseif iszero(w.b) && w.a == a #L_6^t
-        A = _BandedMatrix(-(a:∞)', ∞, 0,0)
+        A = Diagonal(-(a:∞))
         ApplyQuasiMatrix(*, JacobiWeight(a-1,w.b) .* Jacobi(a-1,b+1), A)
     elseif w.a == a && w.b == b # L_1^t
-        A = _BandedMatrix((-2*(1:∞))', ∞, 1,-1)
+        A = _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1)
         ApplyQuasiMatrix(*, JacobiWeight(a-1,b-1) .* Jacobi(a-1, b-1), A)
     elseif iszero(w.a)
         W = (JacobiWeight(w.a, b-1) .* Jacobi(a+1, b-1)) \ (D * (JacobiWeight(w.a,b) .* S))
@@ -373,13 +369,13 @@ function \(L::Legendre, WS::WeightedBasis{Bool,JacobiWeight{Bool},Jacobi{Bool}})
     w,S = WS.args
     if w.b && w.a
         @assert S.b && S.a
-        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))', Zeros(1,∞), (-(2:2:∞)./(3:2:∞))'), ∞, 2,0)
+        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))', Zeros(1,∞), (-(2:2:∞)./(3:2:∞))'), ℵ₀, 2,0)
     elseif w.b && !w.a
         @assert S.b && !S.a
-        _BandedMatrix(Ones{eltype(L)}(2,∞), ∞, 1,0)
+        Bidiagonal(Ones{eltype(L)}(∞), Ones{eltype(L)}(∞), :L)
     elseif !w.b && w.a
         @assert !S.b && S.a
-        _BandedMatrix(Vcat(Ones{eltype(L)}(1,∞),-Ones{eltype(L)}(1,∞)), ∞, 1,0)
+        Bidiagonal(Ones{eltype(L)}(∞), -Ones{eltype(L)}(1,∞), :L)
     else
         error("Not implemented")
     end

src/lanczos.jl

@@ -177,8 +177,8 @@ end
 
 getindex(A::LanczosRecurrence, I::Integer) = _lanczos_getindex(A, I)
 getindex(A::LanczosRecurrence, I::AbstractVector) = _lanczos_getindex(A, I)
-getindex(K::LanczosRecurrence, k::AbstractInfUnitRange) = view(K, k)
-getindex(K::SubArray{<:Any,1,<:LanczosRecurrence}, k::AbstractInfUnitRange) = view(K, k)
+getindex(K::LanczosRecurrence, k::InfRanges{<:Integer}) = view(K, k)
+getindex(K::SubArray{<:Any,1,<:LanczosRecurrence}, k::InfRanges{<:Integer}) = view(K, k)
 
 
 struct LanczosPolynomial{T,XX,WW,Weight,Basis} <: OrthogonalPolynomial{T}

src/normalized.jl

@@ -12,7 +12,7 @@ end
 
 NormalizationConstant(P::AbstractQuasiMatrix{T}) where T = NormalizationConstant{T,typeof(P)}(P)
 
-size(K::NormalizationConstant) = (∞,)
+size(K::NormalizationConstant) = (ℵ₀,)
 
 # How we populate the data
 # function _normalizationconstant_fill_data!(K::NormalizationConstant, J::Union{BandedMatrix,Symmetric{<:Any,BandedMatrix},Tridiagonal,SymTridiagonal}, inds)
@@ -77,9 +77,9 @@ _p0(Q::Normalized) = Q.scaling[1]
 # q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
 
 function recurrencecoefficients(Q::Normalized)
-    A,B,C = recurrencecoefficients(Q.P)
-    h = Q.scaling
-    h[2:∞] ./ h .* A, h[2:∞] ./ h .* B, Vcat(zero(eltype(Q)), h[3:∞] ./ h .* C[2:∞])
+    X = jacobimatrix(Q.P)
+    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
+    inv.(sqrt.(b .* c)), -(a ./ sqrt.(b .* c)), Vcat(zero(eltype(Q)), sqrt.(b .* c) ./ sqrt.(b[2:end] .* c[2:end]))
 end
 
 # x * p[n] = c[n-1] * p[n-1] + a[n] * p[n] + b[n] * p[n+1]
@@ -88,12 +88,12 @@ end
 
 # q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
 # q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
-function jacobimatrix(Q::Normalized)
-    X = jacobimatrix(Q.P)
-    a,b = X[band(0)], X[band(-1)]
-    h = Q.scaling
-    Symmetric(_BandedMatrix(Vcat(a', (b .* h ./ h[2:end])'), ∞, 1, 0), :L)
+
+function symtridagonalize(X)
+    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
+    SymTridiagonal(a, sqrt.(b .* c))
 end
+jacobimatrix(Q::Normalized) = symtridagonalize(jacobimatrix(Q.P))
 
 orthogonalityweight(Q::Normalized) = orthogonalityweight(Q.P)
 singularities(Q::Normalized) = singularities(Q.P)
@@ -131,7 +131,7 @@ _mul_arguments(Q::QuasiAdjoint{<:Any,<:Normalized}) = arguments(ApplyLayout{type
 
 # table stable identity if A.P == B.P
 @inline _normalized_ldiv(An, C, Bn) = An \ (C * Bn)
-@inline _normalized_ldiv(An, C::Eye{T}, Bn) where T = FillArrays.SquareEye{promote_type(eltype(An),T,eltype(Bn))}(∞)
+@inline _normalized_ldiv(An, C::Eye{T}, Bn) where T = FillArrays.SquareEye{promote_type(eltype(An),T,eltype(Bn))}(ℵ₀)
 @inline copy(L::Ldiv{<:NormalizedBasisLayout,<:NormalizedBasisLayout}) = _normalized_ldiv(Diagonal(L.A.scaling), L.A.P \ L.B.P, Diagonal(L.B.scaling))
 @inline copy(L::Ldiv{Lay,<:NormalizedBasisLayout}) where Lay = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A, L.B))
 @inline copy(L::Ldiv{<:NormalizedBasisLayout,Lay}) where Lay = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A, L.B))

src/stieltjes.jl

@@ -14,12 +14,12 @@ end
 
 @simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wT::WeightedBasis{<:Any,<:ChebyshevTWeight,<:ChebyshevT}) 
     T = promote_type(eltype(H), eltype(wT))
-    ApplyQuasiArray(*, ChebyshevU{T}(), _BandedMatrix(Fill(-convert(T,π),1,∞), ∞, -1, 1))
+    ApplyQuasiArray(*, ChebyshevU{T}(), _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1))
 end
 
 @simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval}, wU::WeightedBasis{<:Any,<:ChebyshevUWeight,<:ChebyshevU}) 
     T = promote_type(eltype(H), eltype(wU))
-    ApplyQuasiArray(*, ChebyshevT{T}(), _BandedMatrix(Fill(convert(T,π),1,∞), ∞, 1, -1))
+    ApplyQuasiArray(*, ChebyshevT{T}(), _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1))
 end
 
 ###