Recurrence coefficients

Author: DanielVandH

src/ClassicalOrthogonalPolynomials.jl

@@ -325,8 +325,7 @@ function gaussradau(P::Monic{T}, n::Integer, endpt) where {T}
     J = jacobimatrix(P.P, n + 1)
     α, β = diagonaldata(J), supdiagonaldata(J)
     endpt = T(endpt)
-    p0 = P[endpt, n] # πₙ₋₁
-    p1 = P[endpt, n+1] # πₙ. Could be faster by computing p1 from p0 and πₙ₋₂, but the cost is tiny relative to eigen()
+    p0, p1 = P[endpt,n:n+1]
     a = (endpt - β[end]^2 * p0 / p1)::T
     α′ = vcat(@view(α[begin:end-1]), a)
     J′ = SymTridiagonal(α′, β)
@@ -347,8 +346,8 @@ function gausslobatto(P::Monic{T}, n::Integer) where {T}
     a, b = axes(P, 1)[begin], axes(P, 1)[end]
     J = jacobimatrix(P.P, n + 2)
     α, β = diagonaldata(J), supdiagonaldata(J)
-    p0a, p0b = P[a, n+1], P[b, n+1]
-    p1a, p1b = P[a, n+2], P[b, n+2]
+    p0a, p0b = P[a, (n+1):(n+2)]
+    p1a, p1b = P[b, (n+1):(n+2)]
     # Solve Eq. 3.1.2.8 
     Δ = p1a * p0b - p1b * p0a # This could underflow/overflow for large n 
     aext = (a * p1a * p0b - b * p1b * p0a) / Δ

src/monic.jl

@@ -1,13 +1,17 @@
 struct Monic{T,OPs<:AbstractQuasiMatrix{T},NL} <: OrthogonalPolynomial{T}
     P::Normalized{T,OPs,NL}
-    # X::AbstractMatrix{T} # Should X be stored? Probably not. 
-    # X would be the Jacobi matrix of the normalised polynomials, not for the monic polynomials.
+    α::AbstractVector{T} # diagonal of Jacobi matrix of P
+    β::AbstractVector{T} # squared supdiagonal of Jacobi matrix of P
 end
-# Will need to figure out what this should be exactly. 
-# Consider this internal for now until it stabilises.
 
 Monic(P::AbstractQuasiMatrix) = Monic(Normalized(P))
-Monic(P::Monic) = P
+function Monic(P::Normalized)
+    X = jacobimatrix(P)
+    α = diagonaldata(X)
+    β = supdiagonaldata(X)
+    return Monic(P, α, β.^2)
+end
+Monic(P::Monic) = Monic(P.P, P.α, P.β)
 
 Normalized(P::Monic) = P.P
 
@@ -20,23 +24,15 @@ _p0(::Monic{T}) where {T} = one(T)
 show(io::IO, P::Monic) = print(io, "Monic($(P.P.P))")
 show(io::IO, ::MIME"text/plain", P::Monic) = show(io, P)
 
-function getindex(P::Monic{T}, x::Number, n::Int)::T where {T}
-    # TODO: Rewrite this to be more efficient using forwardrecurrence!
-    p0 = _p0(P)
-    n == 1 && return p0
-    t = convert(T, x)
-    J = jacobimatrix(P.P, n)
-    α = diagonaldata(J)
-    β = supdiagonaldata(J)
-    p1 = ((t - α[1]) * p0)::T
-    n == 2 && return p1
-    for i in 2:(n-1)
-        _p1 = p0::T
-        p0 = p1::T
-        p1 = ((t - α[i]) * p1 - β[i-1]^2 * _p1)::T
-    end
-    return p1
+function recurrencecoefficients(P::Monic{T}) where {T}
+    α = P.α
+    β = P.β
+    return _monicrecurrencecoefficients(α, β) # function barrier 
+end
+function _monicrecurrencecoefficients(α::AbstractVector{T}, β) where {T}
+    A = Ones{T}(∞)
+    B = -α
+    C = Vcat(zero(T), β)
+    return A, B, C
 end
-# Should a method be written that makes this more efficient when requesting multiple n?
-
 

test/test_monic.jl

@@ -1,6 +1,6 @@
 using ClassicalOrthogonalPolynomials
 using Test
-using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight
+using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight, recurrencecoefficients
 
 @testset "Basic definition" begin
     P1 = Legendre()
@@ -15,7 +15,7 @@ using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight
     @test sprint(show, MIME"text/plain"(), P3) == "Monic(Legendre())"
 end
 
-@testset "getindex" begin
+@testset "evaluation" begin
     function _pochhammer(x, n)
         y = one(x)
         for i in 0:(n-1)
@@ -53,6 +53,6 @@ end
         @test Q[0.17, 3] ≈ _P(0.17, 2)
         @test Q[0.4, 17] ≈ _P(0.4, 16)
         @test Q[0.9, 21] ≈ _P(0.9, 20)
-        @inferred Q[0.4, 5]
+        # @inferred Q[0.2, 5] # no longer inferred
     end
 end