Support more Lanczos singularitities (#7)

* Support more Lanczos singularitities

* v0.1.2

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/jacobi.jl

@@ -68,6 +68,7 @@ for op in (:+, :-, :*)
         singularitiesbroadcast(::typeof($op), ::NoSingularities, L::LegendreWeight) = L
     end
 end
+singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
 singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
 
 _parent(::NoSingularities) = NoSingularities()
@@ -78,8 +79,8 @@ singularitiesbroadcast(F::Function, G::Function, V::SubQuasiArray, K) = singular
 singularitiesbroadcast(F, V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(F, map(_parent,V)...)[_parentindices(V...)...]
 
 
-singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::JacobiWeight) = b
-singularitiesbroadcast(::typeof(*), a::JacobiWeight, ::LegendreWeight) = a
+singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b
+singularitiesbroadcast(::typeof(*), a::AbstractJacobiWeight, ::LegendreWeight) = a
 
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end
 

test/test_lanczos.jl

@@ -183,4 +183,19 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout
         @test (pϕ.P' * (ϕw .* pϕ.P))[1:3,1:3] ≈ [2 -0.5 0; -0.5 2 -0.5; 0 -0.5 2]
         @test (pϕ' * (ϕw .* pϕ))[1:5,1:5] ≈ I
     end
+
+    @testset "weighted Chebyshev" begin
+        T = ChebyshevT()
+        x = axes(T,1)
+        t = 2
+        Q = LanczosPolynomial((t .- x).^(-1) .* ChebyshevWeight())
+        # perturbation of Toeplitz
+        @test jacobimatrix(Q)[3:10,3:10] ≈ Symmetric(BandedMatrix(1 => Fill(0.5,7)))
+
+        Q = LanczosPolynomial((t .- x) .^ (-1))
+        @test Q[0.1,2] ≈ -0.1327124082839674
+
+        R = LanczosPolynomial((t .- x).^(-1/2) .* ChebyshevWeight())
+        @test R[0.1,2] ≈ -0.03269577983003056
+    end
 end