juliaapproximation / classicalorthogonalpolynomials.jl Goto Github PK

View Code? Open in Web Editor NEW

38.0 38.0 6.0 1.05 MB

A Julia package for classical orthogonal polynomials and expansions

License: MIT License

Julia 90.65% Mathematica 9.35%

classicalorthogonalpolynomials.jl's People

Contributors

Stargazers

Watchers

Forkers

marcofasondini putianyi889 ioannispapapadopoulos jamiecjx tsgut danielvandh

classicalorthogonalpolynomials.jl's Issues

no method matching combine_axes() when indexing a transposed \(A, B)

I doubt this is specifically a ClassicalOrthogonalPolynomials.jl issue but I can't seem to find what types are actually leading to the error because my debugger keeps crashing trying to step through. e.g. similar code with just some LazyBandedMatrices works fine...

julia> using ClassicalOrthogonalPolynomials

julia> A = Jacobi(1.0, 1.0); B = Jacobi(0.0, 1.0);

julia> (A \ B)[1:10, 1:10]; # fine

julia> (A \ B)'[1:10, 1:10]; # fine
10×10 BandedMatrix{Float64} with bandwidths (1, 0):
  1.0    ⋅      ⋅          ⋅       ⋅     ⋅          ⋅          ⋅        ⋅         ⋅
 -0.5   0.75    ⋅          ⋅       ⋅     ⋅          ⋅          ⋅        ⋅         ⋅
   ⋅   -0.5    0.666667    ⋅       ⋅     ⋅          ⋅          ⋅        ⋅         ⋅
   ⋅     ⋅    -0.5        0.625    ⋅     ⋅          ⋅          ⋅        ⋅         ⋅
   ⋅     ⋅      ⋅        -0.5     0.6    ⋅          ⋅          ⋅        ⋅         ⋅
   ⋅     ⋅      ⋅          ⋅     -0.5   0.583333    ⋅          ⋅        ⋅         ⋅
   ⋅     ⋅      ⋅          ⋅       ⋅   -0.5        0.571429    ⋅        ⋅         ⋅
   ⋅     ⋅      ⋅          ⋅       ⋅     ⋅        -0.5        0.5625    ⋅         ⋅
   ⋅     ⋅      ⋅          ⋅       ⋅     ⋅          ⋅        -0.5      0.555556   ⋅
   ⋅     ⋅      ⋅          ⋅       ⋅     ⋅          ⋅          ⋅      -0.5       0.55

julia> (Normalized(A) \ B)'[1:10, 1:10]
ERROR: MethodError: no method matching combine_axes()

Closest candidates are:
  combine_axes(::Any)
   @ Base broadcast.jl:525
  combine_axes(::Any, ::Any)
   @ Base broadcast.jl:524
  combine_axes(::Any, ::Any...)
   @ Base broadcast.jl:523

Stacktrace:
  [1] _axes
    @ .\broadcast.jl:236 [inlined]
  [2] axes
    @ .\broadcast.jl:234 [inlined]
  [3] copyto!
    @ .\broadcast.jl:992 [inlined]
  [4] copyto!
    @ .\broadcast.jl:967 [inlined]
  [5] _BandedMatrix(::LazyBandedMatrices.BroadcastBandedLayout{…}, V::SubArray{…})
    @ LazyBandedMatrices C:\Users\User\.julia\packages\LazyBandedMatrices\59umZ\src\LazyBandedMatrices.jl:575
  [6] BandedMatrices.BandedMatrix(A::SubArray{Float64, 2, LinearAlgebra.Adjoint{…}, Tuple{…}, false})
    @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\dec3g\src\banded\BandedMatrix.jl:245
  [7] sub_materialize(::LazyBandedMatrices.BroadcastBandedLayout{…}, V::SubArray{…}, ::Tuple{…})
    @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\dec3g\src\generic\indexing.jl:28
  [8] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:127 [inlined]
  [9] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:128 [inlined]
 [10] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:134 [inlined]
 [11] getindex(A::LinearAlgebra.Adjoint{Float64, LazyArrays.BroadcastMatrix{…}}, kr::UnitRange{Int64}, jr::UnitRange{Int64})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:149
 [12] top-level scope
    @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> (Normalized(A) \ B)[1:10, 1:10] # but this is fine
10×10 BandedMatrix{Float64} with bandwidths (0, 1):
 1.1547  -0.57735     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅       0.774597  -0.516398    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅         0.617213  -0.46291     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅          ⋅         0.527046  -0.421637    ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅          ⋅          ⋅         0.467099  -0.389249    ⋅          ⋅          ⋅          ⋅
  ⋅        ⋅          ⋅          ⋅          ⋅         0.423659  -0.363137    ⋅          ⋅          ⋅
  ⋅        ⋅          ⋅          ⋅          ⋅          ⋅         0.39036   -0.341565    ⋅          ⋅
  ⋅        ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         0.363803  -0.323381    ⋅
  ⋅        ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         0.341993  -0.307794
  ⋅        ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         0.323669

julia> typeof(Normalized(A) \ B)
BroadcastMatrix{Float64, typeof(\), Tuple{Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}, Bidiagonal{Float64, BroadcastVector{Float64, typeof(/), Tuple{InfStepRange{Float64, Float64}, InfStepRange{Float64, Float64}}}, BroadcastVector{Float64, typeof(/), Tuple{InfStepRange{Float64, Float64}, InfStepRange{Float64, Float64}}}}}} (alias for LazyArrays.BroadcastArray{Float64, 2, typeof(\), Tuple{LazyArrays.Accumulate{Float64, 1, typeof(*), Array{Float64, 1}, AbstractArray{Float64, 1}}, LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastArray{Float64, 1, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastArray{Float64, 1, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}})

julia> typeof((Normalized(A) \ B)')
LinearAlgebra.Adjoint{Float64, LazyArrays.BroadcastMatrix{Float64, typeof(\), Tuple{LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}, LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}

How to get BigFloat accurate conversion operators?

For OPs on y⁴ = 1 - x⁴ I need conversion operators between ChebyshevT() and certain semiclassical OPs, e.g.,

T = ChebyshevT{BigFloat}()
x = axes(T,1)
P₁ = LanczosPolynomial( (1 .+ x.^2).^(0.5) )
R₁ = P₁\T

This is accurate to only Float64 precision:

xv = big"0.1"
T[xv,3] - ( P₁[xv,1]*R₁[1,3] + P₁[xv,3]*R₁[3,3] )
-6.05859383536196727705447569248701020544252033186250536144991618812281214296214e-17

Another example:

P = Legendre{BigFloat}()
R = P\T

The entry R[3,3] is 4/3 but it's accurate to only Float64 precision:

R[3,3]
1.333333333333333250903485233170513237354072211061085640237818614115573313606597

NaN results from Jacobi expansion

julia> @time Chebyshev()\JacobiWeight(0,0.5)
  0.027141 seconds (3.08 k allocations: 3.701 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
  0.9003163161786569
  0.600210877394969
 -0.12004217544451247
  0.05144664659444724
 -0.028581470311092174
  ⋮

julia> @time Jacobi(-0.5,-0.5)\JacobiWeight(0,0.5)
 14.685838 seconds (3.11 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
 NaN
 NaN
 NaN
 NaN
 NaN
   ⋮

julia> @time Legendre()\JacobiWeight(0,0.5)
 15.032195 seconds (3.09 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
  0.9428090415820628
  0.5656854249492391
 -0.13468700594029615
  0.0628539361054728
 -0.036732819801901045
  ⋮

julia> @time Jacobi(0,0)\JacobiWeight(0,0.5)
 14.450656 seconds (3.11 k allocations: 4.702 MiB)
ℵ₀-element LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
 NaN
 NaN
 NaN
 NaN
 NaN
   ⋮

Can't print slices of orthogonal polynomials

While the following variables can be defined, printing them to REPL errors:

using ClassicalOrthogonalPolynomials
Chebyshev()[0.3, :] # errors when printing
Jacobi(0.0, 0.0)[0.3, :] # errors when printing
Laguerre(0)[0.3, :] # errors when printing
# etc

with errors like

julia> Chebyshev()[0.3, :]
ℵ₀-element view(::ChebyshevT{Float64}, 0.3, :) with eltype Float64 with indices OneToInf():
ERROR: MethodError: no method matching isassigned(::ChebyshevT{Float64}, ::Float64, ::Int64)

Closest candidates are:
  isassigned(::Array, ::Int64...)
   @ Base array.jl:266
  isassigned(::LinearAlgebra.UnitLowerTriangular, ::Int64, ::Int64)
   @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:226
  isassigned(::Base.ReshapedArray{T, N}, ::Int64...) where {T, N}
   @ Base reshapedarray.jl:235
  ...

Stacktrace:
  [1] isassigned
    @ .\subarray.jl:362 [inlined]
  [2] isassigned(::SubArray{Float64, 1, ChebyshevT{Float64}, Tuple{Float64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, ::Int64, ::Int64)
    @ Base .\multidimensional.jl:1575
  [3] alignment(io::IOContext{Base.TTY}, X::AbstractVecOrMat, rows::Vector{Int64}, cols::Vector{Int64}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Int64)
    @ Base .\arrayshow.jl:68
  [4] _print_matrix(io::IOContext{Base.TTY}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{Int64}, colsA::UnitRange{Int64})
    @ Base .\arrayshow.jl:207
  [5] print_matrix(io::IOContext{Base.TTY}, X::SubArray{Float64, 1, ChebyshevT{Float64}, Tuple{Float64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base .\arrayshow.jl:171
  [6] print_matrix
    @ .\arrayshow.jl:171 [inlined]
  [7] print_array
    @ .\arrayshow.jl:358 [inlined]
  [8] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::SubArray{Float64, 1, ChebyshevT{Float64}, Tuple{Float64, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false})
    @ Base .\arrayshow.jl:399
  [9] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:273
 [10] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [11] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:259
 [12] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:278
 [13] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [14] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [15] invokelatest
    @ .\essentials.jl:889 [inlined]
 [16] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:237
 [17] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [18] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [19] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [20] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:150
 [21] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:515
 [22] with_logger
    @ .\logging.jl:627 [inlined]
 [23] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [24] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [25] invokelatest(::Any)
    @ Base .\essentials.jl:889
 [26] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:34

Judging from the error, would it be reasonable to define

using ClassicalOrthogonalPolynomials: OrthogonalPolynomial
Base.isassigned(P::OrthogonalPolynomial, x, n) = (x ∈ axes(P, 1)) && (n ∈ axes(P, 2))

The printing works with this:

julia> Chebyshev()[0.3, :]' # transposing to reduce vertical space
1×ℵ₀ adjoint(view(::ChebyshevT{Float64}, 0.3, :)) with eltype Float64 with indices Base.OneTo(1)×OneToInf():
 1.0  0.3  -0.82  -0.792  0.3448  0.99888  0.254528  -0.846163  -0.762226  0.388828  …

julia> Jacobi(0.0, 0.0)[0.3, :]'
1×ℵ₀ adjoint(view(::Jacobi{Float64}, 0.3, :)) with eltype Float64 with indices Base.OneTo(1)×OneToInf():
 1.0  0.3  -0.365  -0.3825  0.0729375  0.345386  0.129181  -0.224073  -0.239075  …

julia> Laguerre(0)[0.3, :]'
1×ℵ₀ adjoint(view(::Laguerre{Float64}, 0.3, :)) with eltype Float64 with indices Base.OneTo(1)×OneToInf():
 1.0  0.7  0.445  0.2305  0.0523375  -0.0933327  -0.210058  -0.301106  -0.369481  …

Can't compute scalar times squared jacobmatrix

Expressions like 1jacobimatrix(P)^2 seem to never complete. Not sure if it's because of something in this package or not - tried getting a specific matrix using just LazyArrays or LazyBandedMatrices but can't seem to get a better MWE of this.

julia> using ClassicalOrthogonalPolynomials

julia> P = Jacobi(0.2, 0.5);

julia> X = jacobimatrix(P);

julia> 1X; # works fine

julia> 2X # works fine
ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}, Int64}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}, Int64}}} with indices OneToInf()×OneToInf():
 0.222222  0.720721    ⋅          ⋅           ⋅          …
 1.48148   0.0330969  0.821202    ⋅           ⋅
  ⋅        1.24209    0.0133376  0.868385     ⋅
  ⋅         ⋅         1.16261    0.00720535  0.895841
  ⋅         ⋅          ⋅         1.12256     0.00451176
  ⋅         ⋅          ⋅          ⋅          1.09837     …
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅          …
  ⋅         ⋅          ⋅          ⋅           ⋅
  ⋅         ⋅          ⋅          ⋅           ⋅
 ⋮                                                       ⋱

julia> 1X^2 # never ends
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [6] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [7] macro expansion
    @ .\broadcast.jl:1004 [inlined]
  [8] macro expansion
    @ .\simdloop.jl:77 [inlined]
  [9] copyto!
    @ .\broadcast.jl:1003 [inlined]
 [10] copyto!
    @ .\broadcast.jl:956 [inlined]
 [11] copy
    @ .\broadcast.jl:928 [inlined]
 [12] materialize
    @ .\broadcast.jl:903 [inlined]
 [13] broadcast(::typeof(*), ::Int64, ::LazyBandedMatrices.Tridiagonal{…})
    @ Base.Broadcast .\broadcast.jl:841
 [14] broadcasted(::LazyArrays.LazyArrayStyle{…}, ::typeof(*), a::Int64, b::LazyArrays.ApplyArray{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\mul.jl:284
 [15] broadcasted
    @ .\broadcast.jl:1347 [inlined]
 [16] broadcast_preserving_zero_d
    @ .\broadcast.jl:891 [inlined]
 [17] *(A::Int64, B::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{…}})
    @ Base .\arraymath.jl:21
 [18] top-level scope
    @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> typeof(X^2)
LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}

julia> typeof(X)
LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Int64, Int64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}

extrapolation syntax?

It's sometimes helpful to use an OP on a domain bigger than the support of its weight. What syntax to use?

We are already using Base.unsafe_getindex for extrapolation at a point, so we could support, e.g.,

T_ext = Base.unsafe_getindex(Chebyshev(), Inclusion(-2..2), :)

though perhaps that's not explicit enough.

There's also the question of what type T_ext should be. SuperQuasiArray (as opposed to SubQuasiArray)?

Cholesky Jacobi matrices are insanely slow

The culprit is that UX is a lazy multiplication so the following line doesn't do anything:

ClassicalOrthogonalPolynomials.jl/src/choleskyQR.jl

Line 104 in fb28b6f

resizedata!(J.UX,inds[end]+1,inds[end]+1)

Here is the slow code:

julia> r = 0.5
0.5

julia> lmin, lmax = (1-sqrt(r))^2,  (1+sqrt(r))^2
(0.08578643762690492, 2.914213562373095)

julia> U = chebyshevu(lmin..lmax);

julia> Q = OrthogonalPolynomial(x -> 1/(2π) * sqrt((lmax-x)*(x-lmin))/(x*r), U);

julia> J = jacobimatrix(Q);

julia> @time J[1:1000,1:1000]
 36.732166 seconds (178.69 k allocations: 1.624 GiB, 0.29% gc time)
1000×1000 BandedMatrix{Float64} with bandwidths (1, 1):
 1.0       0.707107   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅   …   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
 0.707107  1.5       0.707107   ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅        0.707107  1.5       0.707107   ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅         ⋅        0.707107  1.5       0.707107   ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅        0.707107  1.5       0.707107   ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        0.707107  1.5       0.707107   ⋅    ⋅    ⋅    ⋅    ⋅   …   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅ 
 ⋮                                                 ⋮                                  ⋮         ⋱                 ⋮                             ⋮                                       
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅   …   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.707107  1.5       0.707107   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅        0.707107  1.5       0.707107   ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅        0.707107  1.5       0.707107   ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅        0.707107  1.5       0.707107
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅       ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅         ⋅         ⋅         ⋅        0.707107  1.5

Unify LegendreWeight and QuasiOnes

They are the same thing though the Weight interface might be an issue

Make adaptive transform a Vcat instead of a cached vector?

Sometimes we use the transform to build operators, and the bandwidth is determined by the number of non-zeros. In this case, cached vectors are not so useful as anytime we access them they grow, and the bandwidth of the operators grows.

But is there any reason we would need coefficient vectors to be mutable (in their zeros) by default? A user could call Base.copymutable(P \ f) if they really wanted to mutate.

Ultraspherical{BigFloat} expansion broken

Solutions:

Use GenericLinearAlgebra.jl
Use Chebyshev->Ultraspherical transform

julia> C = Ultraspherical{BigFloat}(2)
Ultraspherical{BigFloat, Int64}

julia> x = axes(C,2)
OneToInf()

julia> C \ exp.(x)^C

julia> x = axes(C,1)
Inclusion(-1.0..1.0 (Chebyshev))

julia> C \ exp.(x)
ERROR: MethodError: no method matching eigen!(::SymTridiagonal{BigFloat, Vector{BigFloat}})
Closest candidates are:
  eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{Float32, Float64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:148
  eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{ComplexF32, ComplexF64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:171
  eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at /Users/sheehanolver/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:427
  ...
Stacktrace:
  [1] eigen(A::SymTridiagonal{BigFloat, Vector{BigFloat}})
    @ LinearAlgebra ~/Projects/julia-1.6/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/tridiag.jl:281
  [2] eigen(A::LazyBandedMatrices.SymTridiagonal{BigFloat, FillArrays.Zeros{BigFloat, 1, Tuple{Base.OneTo{Int64}}}, Vector{BigFloat}})
    @ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/IX6Zy/src/tridiag.jl:723
  [3] golubwelsch(X::LazyBandedMatrices.Tridiagonal{BigFloat, Vector{BigFloat}, FillArrays.Zeros{BigFloat, 1, Tuple{Base.OneTo{Int64}}}, Vector{BigFloat}})
    @ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:287
  [4] golubwelsch(V::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false})
    @ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:292
  [5] factorize(L::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false})
    @ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:304
  [6] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64})
    @ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:255
  [7] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
    @ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:256
  [8] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, QuasiArrays.SubQuasiArray{BigFloat, 2, Ultraspherical{BigFloat, Int64}, Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}, Base.OneTo{Int64}}, false}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
    @ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:262
  [9] materialize
    @ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:22 [inlined]
 [10] ldiv
    @ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:86 [inlined]
 [11] \
    @ ~/Projects/QuasiArrays.jl/src/matmul.jl:34 [inlined]
 [12] adaptivetransform_ldiv(A::Ultraspherical{BigFloat, Int64}, f::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
    @ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:92
 [13] transform_ldiv(A::Ultraspherical{BigFloat, Int64}, f::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Infinities.InfiniteCardinal{0}})
    @ ClassicalOrthogonalPolynomials ~/Projects/ClassicalOrthogonalPolynomials.jl/src/ClassicalOrthogonalPolynomials.jl:62
 [14] transform_ldiv(A::Ultraspherical{BigFloat, Int64}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
    @ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:256
 [15] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.BasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, Ultraspherical{BigFloat, Int64}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
    @ ContinuumArrays ~/Projects/ContinuumArrays.jl/src/bases/bases.jl:262
 [16] materialize(M::ArrayLayouts.Ldiv{ContinuumArrays.BasisLayout, LazyArrays.BroadcastLayout{typeof(exp)}, Ultraspherical{BigFloat, Int64}, QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:22
 [17] ldiv
    @ ~/.julia/packages/ArrayLayouts/3lnt1/src/ldiv.jl:86 [inlined]
 [18] \(A::Ultraspherical{BigFloat, Int64}, B::QuasiArrays.BroadcastQuasiVector{BigFloat, typeof(exp), Tuple{Inclusion{BigFloat, ChebyshevInterval{BigFloat}}}})
    @ QuasiArrays ~/Projects/QuasiArrays.jl/src/matmul.jl:34
 [19] top-level scope
    @ REPL[130]:1

Shifted Legendre{BigFloat}() error

I'm seeing some strange behavior for Legendre polynomials with BigFloat precision shifted to 0..1. Quite possible that something about my syntax isn't how it is intended to be used but it might also be bug.

Here's a simple example that reproduces the error:

julia> # float64

julia> P = Legendre()[affine(Inclusion(0..1), axes(Legendre(),1)), :]
Legendre{Float64} affine mapped to 0..1

julia> x = axes(P,1)
Inclusion(0..1)

julia> P \ sin.(x)
ℵ₀-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
  0.4596976941318605
  0.42791899124296007
 -0.03924363301542014
 -0.007212191228055754
  0.0002821450243386184
  2.8742703093836133e-5
 -7.146539529763204e-7
 -5.0363463994812964e-8
  9.178337836736593e-10
  4.944523575181289e-11
 -7.110265998998354e-13
 -3.106481507574827e-14
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  ⋮

julia> # bigfloat

julia> Pbig = Legendre{BigFloat}()[affine(Inclusion(BigInt(0)..BigInt(1)), axes(Legendre{BigFloat}(),1)), :]
Legendre{BigFloat} affine mapped to 0..1

julia> xbig = axes(Pbig,1)
Inclusion(0..1)

julia> Pbig \ sin.(xbig)
ERROR: cannot resize array with shared data
Stacktrace:
 [1] _deleteend! at ./array.jl:901 [inlined]
 [2] resize! at ./array.jl:1090 [inlined]
 [3] chop!(::Array{BigFloat,1}, ::BigFloat) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:71
 [4] adaptivetransform_ldiv(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:104
 [5] transform_ldiv at /home/timon/.julia/packages/ClassicalOrthogonalPolynomials/YpJSc/src/ClassicalOrthogonalPolynomials.jl:64 [inlined]
 [6] transform_ldiv at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:178 [inlined]
 [7] copy at /home/timon/.julia/packages/ContinuumArrays/D2vwQ/src/bases/bases.jl:180 [inlined]
 [8] materialize(::Ldiv{ContinuumArrays.MappedBasisLayout,LazyArrays.BroadcastLayout{typeof(sin)},SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false},BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}}) at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:22
 [9] ldiv at /home/timon/.julia/packages/ArrayLayouts/7LZ91/src/ldiv.jl:86 [inlined]
 [10] \(::SubQuasiArray{BigFloat,2,Legendre{BigFloat},Tuple{ContinuumArrays.AffineMap{BigFloat,Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}},Inclusion{BigFloat,ChebyshevInterval{BigFloat}}},Base.Slice{InfiniteArrays.OneToInf{Int64}}},false}, ::BroadcastQuasiArray{BigFloat,1,typeof(sin),Tuple{Inclusion{BigFloat,IntervalSets.Interval{:closed,:closed,BigInt}}}}) at /home/timon/.julia/packages/QuasiArrays/Kytgm/src/matmul.jl:34
 [11] top-level scope at REPL[77]:1

julia> # but this appears to work

julia> Pbig[:,1:20] \ sin.(xbig)
20-element Array{BigFloat,1}:
  0.459697694131860282599063392557024291715022853422108373904832662310088291651591
  0.4279189912429598877122041074528643054476364481970978558557517258602675648095019
 -0.03924363301542034337341694172731222575488871592208301467640058895035480424366479
 -0.007212191228055742704936278556706084716636823072871192246094339888763058783659867
  0.000282145024338535280236421222451671910237190823901368171354926761601708489205559
  2.874270309358435249957361731260692515947881520105291315102420408358148420958176e-05
 -7.146539530749001084413538229440980152040709571024099798175596892344455384798109e-07
 -5.036346369111653274248451118754811843127054009898515258792902259125311582439218e-08
  9.178336969399690650953052093298209802341083885853023198941221702769575165208465e-10
  4.944521184304420861236925934788833915496435733298095933603120844076503786764353e-11
 -7.112115015695451912101596871391910849417223830934967825666746546919715622639409e-13
 -3.101024774116237937183881067754829492183211029272491093718995995982405455400059e-14
  3.684185721261002718333096225757674226208193004456196754596387429774500100182845e-16
  1.349202245083837855737893815880041795386546710643032360063629130209568728057782e-17
 -1.365328931605316752701716780091452534277698239440733900315583084467897549345079e-19
 -4.310049800103575884451734511728530188560316941533188009045818391107703584309708e-21
  3.798549690418122116191713913664570774313410824956465599814127761088843678732203e-23
  1.053718404267767678750397050992597701028600476518485629180358929672362729284939e-24
 -8.224287969317306291875699136774388048268971744012539584301297942139600288084304e-27
 -2.035023861243751730600511331372270370404905435888710306620612749522021466923814e-28

Clenshaw errs with degree-0 polynomial

julia> using ClassicalOrthogonalPolynomials

julia> a = (x,α) -> α*(1-x^2)
#3 (generic function with 1 method)

julia> T = ChebyshevT()
ChebyshevT()

julia> x = axes(T, 1)
Inclusion(-1.0 .. 1.0 (Chebyshev))

julia> A = T \ (a.(x, 0.1) .* T) # ✓
ℵ₀×ℵ₀ Clenshaw{Float64} with 3 degree polynomial:
  0.05   0.0    -0.025    ⋅       ⋅     …  
  0.0    0.025   0.0    -0.025    ⋅        
 -0.05   0.0     0.05    0.0    -0.025     
   ⋅    -0.025   0.0     0.05    0.0       
   ⋅      ⋅     -0.025   0.0     0.05      
   ⋅      ⋅       ⋅     -0.025   0.0    …  
   ⋅      ⋅       ⋅       ⋅     -0.025     
   ⋅      ⋅       ⋅       ⋅       ⋅        
   ⋅      ⋅       ⋅       ⋅       ⋅        
   ⋅      ⋅       ⋅       ⋅       ⋅        
  ⋮                                     ⋱  

julia> A = T \ (a.(x, 0.) .* T) # ×
ℵ₀×ℵ₀ Clenshaw{Float64} with 0 degree polynomial:
Error showing value of type ClassicalOrthogonalPolynomials.Clenshaw{Float64, Vector{Float64}, LazyArrays.ApplyArray{Int64, 1, typeof(vcat), Tuple{Int64, FillArrays.Fill{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, FillArrays.Zeros{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Ones{Int64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, FillArrays.Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}:
ERROR: MethodError: no method matching getindex(::Float64, ::UnitRange{Int64})

Closest candidates are:
  getindex(::Number, ::BlockArrays.Block{0})
   @ BlockArrays ~/.julia/packages/BlockArrays/g5q5u/src/blockarrayinterface.jl:2
  getindex(::Number, ::Integer)
   @ Base number.jl:96
  getindex(::Number)
   @ Base number.jl:95
  ...

Stacktrace:
  [1] getindex(M::ClassicalOrthogonalPolynomials.Clenshaw{…}, kr::UnitRange{…}, j::Int64)
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/nNR0P/src/clenshaw.jl:304
  [2] getindex
    @ ~/.julia/packages/ClassicalOrthogonalPolynomials/nNR0P/src/clenshaw.jl:307 [inlined]
  [3] isassigned(::ClassicalOrthogonalPolynomials.Clenshaw{…}, ::Int64, ::Int64)
    @ Base ./multidimensional.jl:1578
  [4] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{…})
    @ Base ./arrayshow.jl:68
  [5] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{…}, colsA::InfiniteArrays.InfUnitRange{…})
    @ Base ./arrayshow.jl:207
  [6] print_matrix(io::IOContext{…}, X::ClassicalOrthogonalPolynomials.Clenshaw{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base ./arrayshow.jl:171
  [7] print_matrix
    @ ./arrayshow.jl:171 [inlined]
  [8] print_array
    @ ./arrayshow.jl:358 [inlined]
  [9] show(io::IOContext{…}, ::MIME{…}, X::ClassicalOrthogonalPolynomials.Clenshaw{…})
    @ Base ./arrayshow.jl:399
 [10] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:273
 [11] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [12] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:259
 [13] display
    @ /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:278 [inlined]
 [14] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:340
 [15] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [16] invokelatest
    @ ./essentials.jl:889 [inlined]
 [17] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{…})
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:315
 [18] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:284
 [19] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [20] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:282
 [21] (::REPL.var"#do_respond#80"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:911
 [22] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [23] invokelatest
    @ ./essentials.jl:889 [inlined]
 [24] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/LineEdit.jl:2656
 [25] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:1312
 [26] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/REPL/src/REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

A \ B never completes when B is a `Normalized` basis

julia> using ClassicalOrthogonalPolynomials

julia> P = Jacobi(2.0, 0.5); Q = Jacobi(3.0, 0.5);

julia> P \ Q; # works fine

julia> Normalized(P) \ Q; # works fine

julia> P \ Normalized(Q); # :(
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [6] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [7] copyto_unaliased!(deststyle::IndexCartesian, dest::LazyArrays.CachedArray{…}, srcstyle::IndexCartesian, src::LinearAlgebra.Diagonal{…})
    @ Base .\abstractarray.jl:1102
  [8] copyto!(dest::LazyArrays.CachedArray{…}, src::LinearAlgebra.Diagonal{…})
    @ Base .\abstractarray.jl:1068
  [9] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:255 [inlined]
 [10] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:259 [inlined]
 [11] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:262 [inlined]
 [12] copyto!(dest::LazyArrays.CachedArray{…}, M::ArrayLayouts.Ldiv{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\triangular.jl:209
 [13] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:21 [inlined]
 [14] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [15] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [16] copy(M::ArrayLayouts.Mul{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:90
 [17] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:127 [inlined]
 [18] mul
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:128 [inlined]
 [19] *(A::LazyArrays.ApplyArray{…}, B::LinearAlgebra.Diagonal{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:219
 [20] _copy_ldiv_mul
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:112 [inlined]
 [21] copy
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:116 [inlined]
 [22] copy
    @ C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:123 [inlined]
 [23] basis_ldiv_size
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:327 [inlined]
 [24] copy
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:323 [inlined]
 [25] copy
    @ C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:125 [inlined]
 [26] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [27] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [28] \(A::Jacobi{Float64}, B::Normalized{Float64, Jacobi{…}, LazyArrays.Accumulate{…}})
    @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
 [29] top-level scope
    @ REPL[27]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> Normalized(P) \ Normalized(Q); # :(
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [6] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [7] copyto_unaliased!(deststyle::IndexCartesian, dest::LazyArrays.CachedArray{…}, srcstyle::IndexCartesian, src::LinearAlgebra.Diagonal{…})
    @ Base .\abstractarray.jl:1102
  [8] copyto!(dest::LazyArrays.CachedArray{…}, src::LinearAlgebra.Diagonal{…})
    @ Base .\abstractarray.jl:1068
  [9] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:255 [inlined]
 [10] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:259 [inlined]
 [11] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:262 [inlined]
 [12] copyto!(dest::LazyArrays.CachedArray{…}, M::ArrayLayouts.Ldiv{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\triangular.jl:209
 [13] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:21 [inlined]
 [14] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [15] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [16] copy(M::ArrayLayouts.Mul{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:90
 [17] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:127 [inlined]
 [18] mul
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:128 [inlined]
 [19] *
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:219 [inlined]
 [20] _normalized_ldiv(An::LinearAlgebra.Diagonal{…}, C::LazyArrays.ApplyArray{…}, Bn::LinearAlgebra.Diagonal{…})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:117
 [21] copy
    @ C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:121 [inlined]
 [22] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [23] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [24] \(A::Normalized{…}, B::Normalized{…})
    @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
 [25] top-level scope
    @ REPL[28]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> Chebyshev() \ Normalized(ChebyshevU()); # :(
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] setindex!
    @ .\subarray.jl:331 [inlined]
  [6] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [7] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [8] copyto_unaliased!(deststyle::IndexCartesian, dest::SubArray{…}, srcstyle::IndexCartesian, src::SubArray{…})
    @ Base .\abstractarray.jl:1102
  [9] copyto!
    @ .\abstractarray.jl:1068 [inlined]
 [10] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:255 [inlined]
 [11] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:259 [inlined]
 [12] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:267 [inlined]
 [13] ldiv!(Y::LazyArrays.CachedArray{…}, A::MatrixFactorizations.QR{…}, B::LinearAlgebra.Diagonal{…})
    @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\factorization.jl:179
 [14] _ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:84 [inlined]
 [15] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
 [16] ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:104 [inlined]
 [17] _ldiv!(dest::LazyArrays.CachedArray{…}, A::BandedMatrices.BandedMatrix{…}, B::LinearAlgebra.Diagonal{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:83
 [18] _ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:83 [inlined]
 [19] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
 [20] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:21 [inlined]
 [21] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [22] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [23] copy(M::ArrayLayouts.Mul{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:90
 [24] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:127 [inlined]
 [25] mul
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:128 [inlined]
 [26] *
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:219 [inlined]
 [27] _copy_ldiv_mul
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:112 [inlined]
 [28] copy
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:116 [inlined]
 [29] copy
    @ C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:123 [inlined]
 [30] basis_ldiv_size
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:327 [inlined]
 [31] copy
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:323 [inlined]
 [32] copy
    @ C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\normalized.jl:125 [inlined]
 [33] materialize(M::ArrayLayouts.Ldiv{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22
 [34] ldiv(A::ChebyshevT{Float64}, B::Normalized{Float64, ChebyshevU{…}, LazyArrays.Accumulate{…}}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98
 [35] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [36] \(A::ChebyshevT{Float64}, B::Normalized{Float64, ChebyshevU{…}, LazyArrays.Accumulate{…}})
    @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
 [37] top-level scope
    @ REPL[29]:1
Some type information was truncated. Use `show(err)` to see complete types.

Make P[0.1,10] allocation free

At the moment we unnecessarily build P[0.1,1:10] but the use of forwardrecurrence would avoid this

Chebyshev() \ exp.(im*x) errors

julia> Chebyshev() \ exp.(im*x)
ERROR: InexactError: Float64(0.5443479390547952 + 0.8388595360647676im)
Stacktrace:
  [1] Real
    @ ./complex.jl:44 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:969 [inlined]
  [4] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{ComplexF64}, soffs::Int64, n::Int64)
    @ Base ./array.jl:250
  [5] unsafe_copyto!
    @ ./array.jl:304 [inlined]
  [6] _copyto_impl!
    @ ./array.jl:327 [inlined]
  [7] copyto!
    @ ./array.jl:314 [inlined]
  [8] copyto!
    @ ./array.jl:339 [inlined]
  [9] circcopy!(dest::Vector{Float64}, src::Vector{ComplexF64})
    @ Base ./multidimensional.jl:1244
 [10] copy1
    @ ~/.julia/packages/AbstractFFTs/0uOAT/src/definitions.jl:50 [inlined]
 [11] *
    @ ~/.julia/packages/AbstractFFTs/0uOAT/src/definitions.jl:220 [inlined]
 [12] \(a::ContinuumArrays.TransformFactorization{Float64, ChebyshevGrid{1, Float64}, FastTransforms.ChebyshevTransformPlan{Float64, 1, Vector{Int32}, false, 1, Int64}}, b::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/2jOj3/src/plans.jl:26
 [13] transform_ldiv
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:257 [inlined]
 [14] transform_ldiv
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:258 [inlined]
 [15] transform_ldiv_if_columns
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:289 [inlined]
 [16] transform_ldiv_if_columns
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:290 [inlined]
 [17] copy
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:291 [inlined]
 [18] #materialize#13
    @ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22 [inlined]
 [19] materialize
    @ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22 [inlined]
 [20] #ldiv#20
    @ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
 [21] ldiv
    @ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
 [22] \
    @ ~/.julia/packages/QuasiArrays/EK26C/src/matmul.jl:34 [inlined]
 [23] adaptivetransform_ldiv(A::ChebyshevT{Float64}, f::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/sveBE/src/adaptivetransform.jl:34
 [24] transform_ldiv
    @ ~/.julia/packages/ClassicalOrthogonalPolynomials/sveBE/src/adaptivetransform.jl:2 [inlined]
 [25] transform_ldiv
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:258 [inlined]
 [26] transform_ldiv_if_columns
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:289 [inlined]
 [27] transform_ldiv_if_columns
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:290 [inlined]
 [28] copy
    @ ~/.julia/packages/ContinuumArrays/2jOj3/src/bases/bases.jl:291 [inlined]
 [29] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.OPLayout, LazyArrays.BroadcastLayout{typeof(exp)}, ChebyshevT{Float64}, QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22
 [30] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.OPLayout, LazyArrays.BroadcastLayout{typeof(exp)}, ChebyshevT{Float64}, QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:22
 [31] ldiv(A::ChebyshevT{Float64}, B::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86
 [32] ldiv
    @ ~/.julia/packages/ArrayLayouts/ph6IB/src/ldiv.jl:86 [inlined]
 [33] \(A::ChebyshevT{Float64}, B::QuasiArrays.BroadcastQuasiVector{ComplexF64, typeof(exp), Tuple{ContinuumArrays.AffineQuasiVector{ComplexF64, ComplexF64, Inclusion{Float64, ChebyshevInterval{Float64}}, ComplexF64}}})
    @ QuasiArrays ~/.julia/packages/QuasiArrays/EK26C/src/matmul.jl:34
 [34] top-level scope
    @ REPL[39]:1

It's because the transform needs to take into account the type of RHS.

Stackoverflow for P'P

Seems to happen specifically for normalized polynomials.

This works:

julia> P = Jacobi(0,1)
Jacobi(0.0, 1.0)

julia> P'P
((inv(ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf())' with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}} with indices OneToInf()×OneToInf()) * (inv(ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
  2.0       -1.0        0.666667  -0.5       0.4      -0.333333   0.285714  -0.25   0.222222  …  
 -1.0        2.0       -1.33333    1.0      -0.8       0.666667  -0.571429   0.5   -0.444444     
  0.666667  -1.33333    2.0       -1.5       1.2      -1.0        0.857143  -0.75   0.666667     
 -0.5        1.0       -1.5        2.0      -1.6       1.33333   -1.14286    1.0   -0.888889     
  0.4       -0.8        1.2       -1.6       2.0      -1.66667    1.42857   -1.25   1.11111      
 -0.333333   0.666667  -1.0        1.33333  -1.66667   2.0       -1.71429    1.5   -1.33333   …  
  0.285714  -0.571429   0.857143  -1.14286   1.42857  -1.71429    2.0       -1.75   1.55556      
 -0.25       0.5       -0.75       1.0      -1.25      1.5       -1.75       2.0   -1.77778      
  ⋮                                                    ⋮                                      ⋱

This fails:

julia> P = Normalized(Jacobi(0,1))
Normalized(Jacobi(0.0, 1.0))

julia> P'P
ERROR: StackOverflowError:
Stacktrace:
     [1] axes
       @ ./abstractarray.jl:74 [inlined]
     [2] axes(D::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/memorylayout.jl:716
     [3] axes
       @ ./abstractarray.jl:74 [inlined]
     [4] _check_mul_axes(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:91
     [5] check_mul_axes
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:93 [inlined]
     [6] instantiate
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:105 [inlined]
     [7] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
     [8] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
     [9] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
    [10] _simplify(#unused#::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293
    [11] simplify(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292
    [12] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
    [13] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:903
--- the last 7 lines are repeated 4041 more times ---
 [28301] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, ArrayLayouts.DiagonalLayout{LazyArrays.LazyLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
 [28302] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
 [28303] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
 [28304] _simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293 [inlined]
 [28305] _tail_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
 [28306] _most_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298 [inlined]
 [28307] _simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296 [inlined]
 [28308] _tail_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
 [28309] _most_simplify(#unused#::Val{false}, args::Tuple{Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298
 [28310] _simplify(::typeof(*), ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
--- the last 3 lines are repeated 1 more time ---
 [28314] simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
 [28315] simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302 [inlined]
 [28316] copy
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:307 [inlined]
 [28317] materialize
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
 [28318] mul
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
 [28319] *(A::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}, B::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/ArrayLayouts.jl:194
 [28320] afoldl(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ Base ./operators.jl:549
 [28321] *(::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ Base ./operators.jl:591
 [28322] _most_simplify(#unused#::Val{true}, args::Tuple{Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:297
 [28323] _simplify(::typeof(*), ::Diagonal{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, ::QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, ::Jacobi{Float64}, ::Vararg{Any})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
 [28324] simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
 [28325] simplify(M::ArrayLayouts.Mul{ContinuumArrays.AdjointBasisLayout{ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}}, ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}, QuasiArrays.QuasiAdjoint{Float64, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
 [28326] copy
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:306 [inlined]
 [28327] materialize
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
 [28328] mul
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
 [28329] *(A::QuasiArrays.QuasiAdjoint{Float64, Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, B::Normalized{Float64, Jacobi{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
       @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:23

A similar but quite possibly unrelated issue also exists for finite sets of polynomials but here it happens even when not normalized:

julia> P = Jacobi(0,1)[:,1:10]
view(Jacobi(0.0, 1.0), Inclusion(-1.0..1.0 (Chebyshev)), 1:10) with eltype Float64

julia> P'P
ERROR: StackOverflowError:
Stacktrace:
     [1] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
     [2] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:899
     [3] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
     [4] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
     [5] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
     [6] _simplify(#unused#::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293
     [7] simplify(::typeof(*), ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292
--- the last 7 lines are repeated 4794 more times ---
 [33566] simplify(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302
 [33567] copy(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ LazyBandedMatrices ~/.julia/packages/LazyBandedMatrices/ZaJ9N/src/LazyBandedMatrices.jl:899
 [33568] materialize(M::ArrayLayouts.Mul{LazyArrays.InvLayout{ArrayLayouts.BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}}, BandedMatrices.BandedColumns{ArrayLayouts.OnesLayout}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109
 [33569] mul(A::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, B::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ ArrayLayouts ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110
 [33570] _twoarg_simplify(#unused#::typeof(*), #unused#::Val{true}, a::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, b::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:295
 [33571] _simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:293 [inlined]
 [33572] _tail_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
 [33573] _most_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298 [inlined]
 [33574] _simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296 [inlined]
 [33575] _tail_simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:299 [inlined]
 [33576] _most_simplify(#unused#::Val{false}, args::Tuple{Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:298
 [33577] _simplify(::typeof(*), ::Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, ::Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, ::LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}, ::Vararg{Any})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:296
--- the last 3 lines are repeated 1 more time ---
 [33581] simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:292 [inlined]
 [33582] simplify
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:302 [inlined]
 [33583] copy
       @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:307 [inlined]
 [33584] materialize
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
 [33585] mul
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
 [33586] *
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:198 [inlined]
 [33587] afoldl(op::typeof(*), a::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, Adjoint{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}, Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{Float64, InfiniteArrays.InfStepRange{Int64, Int64}}}}, LazyArrays.ApplyArray{Float64, 2, typeof(inv), Tuple{LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Float64, Float64}}}}}}}}, bs::BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}})
       @ Base ./operators.jl:548
 [33588] *
       @ ./operators.jl:591 [inlined]
 [33589] _default_materialize(A::LazyArrays.Applied{LazyArrays.DefaultArrayApplyStyle, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}}, #unused#::LazyArrays.Applied{LazyArrays.DefaultArrayApplyStyle, typeof(*), Tuple{Adjoint{Int64, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}, QuasiArrays.QuasiAdjoint{Float64, Jacobi{Float64}}, Jacobi{Float64}, BandedMatrices.BandedMatrix{Int64, Ones{Int64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.OneToInf{Int64}}}})
       @ LazyArrays ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:87
 [33590] copy
       @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:90 [inlined]
 [33591] materialize
       @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:82 [inlined]
 [33592] apply
       @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:79 [inlined]
 [33593] copy(M::ArrayLayouts.Mul{ContinuumArrays.AdjointBasisLayout{ContinuumArrays.SubBasisLayout}, ContinuumArrays.SubBasisLayout, QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}})
       @ ContinuumArrays ~/.julia/packages/ContinuumArrays/sZfkh/src/bases/bases.jl:632
 [33594] materialize
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:109 [inlined]
 [33595] mul
       @ ~/.julia/packages/ArrayLayouts/CNneg/src/mul.jl:110 [inlined]
 [33596] *(A::QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}, B::QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false})
       @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:23

What to do with axes of associated?

Sometimes the Associated OPs contain discrete points not present in the original OPs...

So what do we do for Hilbert transform? It's defined in terms of Associated OPs but we wouldn't want the support to increase.

Perhaps we should properly support restrictions? This would involve introducing QuasiSlice instead of assuming Inclusion is playing that role

Identity mapping wT -> wU

More of a question than an issue. The following works just fine:

Weighted(ChebyshevT())\Weighted(ChebyshevU())

but

Weighted(ChebyshevU())\Weighted(ChebyshevT())

throws an ERROR: OutOfMemoryError() message. Is this expected behaviour? Or do we want something similar to T\U?

Make AbstractOPLayout <: AbstractPolynomialLayout

QuasiArrays.jl now has PolynomialLayout, which currently only supports Inclusion but is designed to represent general polynomials. This package has AbstractOPLayout which should be a special class of polynomial layout.

This helps make sense of how WeightedBasisLayout should work. E.g. equality can be made precise (if w1 and w2 are not polynomial and P and Q are non-zero polynomial quasimatrices then w1 .* P == w2 .* Q iff w1 == w2 and P == Q)

Make p-FEM more versatile

This should work:

julia> P = JacobiWeight(1,1) .* Jacobi(1,1)[:,1:n]
(1-x)^1 * (1+x)^1 on -1..1 .* view(Jacobi(1.0, 1.0), Inclusion(-1.0..1.0 (Chebyshev)), 1:10) with eltype Float64

julia> P'P
(QuasiArrays.QuasiAdjoint{Float64, QuasiArrays.BroadcastQuasiMatrix{Float64, typeof(*), Tuple{JacobiWeight{Int64}, QuasiArrays.SubQuasiArray{Float64, 2, Jacobi{Float64}, Tuple{Inclusion{Float64, ChebyshevInterval{Float64}}, UnitRange{Int64}}, false}}}}) * ((1-x)^1 * (1+x)^1 on -1..1 .* Jacobi(1.0, 1.0)) * (∞×10 BandedMatrix{Int64} with bandwidths (0, 0) with data 1×10 Ones{Int64} with indices OneToInf()×Base.OneTo(10))

Instead we have to do:

julia> P = JacobiWeight(1,1) .* Jacobi(1,1)
(1-x)^1 * (1+x)^1 on -1..1 .* Jacobi(1.0, 1.0)

julia> (P'P)[1:n,1:n]
10×10 BandedMatrix{Float64} with bandwidths (2, 2):
  1.06667    0.0       -0.228571       ⋅             ⋅          ⋅             ⋅             ⋅             ⋅             ⋅ 
  0.0        0.609524   5.55112e-17  -0.203175       ⋅          ⋅             ⋅             ⋅             ⋅             ⋅ 
 -0.228571   0.0        0.457143      0.0          -0.17316     ⋅             ⋅             ⋅             ⋅             ⋅ 
   ⋅        -0.203175   0.0           0.369408      0.0       -0.149184       ⋅             ⋅             ⋅             ⋅ 
   ⋅          ⋅        -0.17316      -1.38778e-17   0.3108     2.77556e-17  -0.130536       ⋅             ⋅             ⋅ 
   ⋅          ⋅          ⋅           -0.149184      0.0        0.268531     -2.77556e-17  -0.115837       ⋅             ⋅ 
   ⋅          ⋅          ⋅             ⋅           -0.130536  -4.16334e-17   0.236501      5.55112e-17  -0.104025       ⋅ 
   ⋅          ⋅          ⋅             ⋅             ⋅        -0.115837      1.38778e-17   0.211352      0.0          -0.0943535
   ⋅          ⋅          ⋅             ⋅             ⋅          ⋅           -0.104025      0.0           0.191066      1.38778e-17
   ⋅          ⋅          ⋅             ⋅             ⋅          ⋅             ⋅           -0.0943535    -1.38778e-17   0.174349
```

Add roots

We can port the code

https://github.com/JuliaApproximation/ApproxFunOrthogonalPolynomials.jl/blob/master/src/roots.jl

I think call it findall(iszero,f).

Conversion operator from LanczosPolynomial

x = Inclusion(-1.0..1)
a = 1.5
ϕ = x.^4 - (a^2 + 1)*x.^2 .+ a^2  
Pϕ = Normalized(LanczosPolynomial(ϕ))
P = Normalized(Legendre())
Cϕ = Pϕ\P

MethodError: copy(::ArrayLayouts.Ldiv{ContinuumArrays.MappedBasisLayout, ContinuumArrays.BasisLayout, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, Legendre{Float64}}) is ambiguous. Candidates:
  copy(P::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:Union{ContinuumArrays.MappedBasisLayout, ContinuumArrays.MappedWeightedBasisLayout}, var"#s12"<:LazyArrays.AbstractLazyLayout, AType, BType}) in ContinuumArrays at C:\Users\mfaso\.julia\packages\ContinuumArrays\gK7Yu\src\bases\bases.jl:88
  copy(P::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:ContinuumArrays.AbstractBasisLayout, var"#s12"<:ContinuumArrays.AbstractBasisLayout, AType, BType}) in ContinuumArrays at C:\Users\mfaso\.julia\packages\ContinuumArrays\gK7Yu\src\bases\bases.jl:71
Possible fix, define
  copy(::ArrayLayouts.Ldiv{var"#s13", var"#s12", AType, BType} where {var"#s13"<:Union{ContinuumArrays.MappedBasisLayout, ContinuumArrays.MappedWeightedBasisLayout}, var"#s12"<:ContinuumArrays.AbstractBasisLayout, AType, BType})

Stacktrace:
  [1] materialize
    @ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22 [inlined]
  [2] ldiv
    @ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
  [3] \
    @ ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34 [inlined]
  [4] _copy_ldiv_ldiv
    @ ~\.julia\packages\LazyArrays\Wldxo\src\linalg\inv.jl:105 [inlined]
  [5] copy
    @ ~\.julia\packages\LazyArrays\Wldxo\src\linalg\inv.jl:107 [inlined]
  [6] copy
    @ ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\normalized.jl:115 [inlined]
  [7] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ContinuumArrays.MappedBasisLayout}, ContinuumArrays.BasisLayout, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Legendre{Float64}})
    @ ArrayLayouts ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22
  [8] ldiv
    @ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
  [9] \
    @ ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34 [inlined]
 [10] \(Q::LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, A::Legendre{Float64})
    @ ClassicalOrthogonalPolynomials ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\lanczos.jl:218
 [11] copy
    @ ~\.julia\packages\ClassicalOrthogonalPolynomials\63Ert\src\normalized.jl:110 [inlined]
 [12] materialize(M::ArrayLayouts.Ldiv{ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ClassicalOrthogonalPolynomials.LanczosLayout}, ClassicalOrthogonalPolynomials.NormalizedBasisLayout{ContinuumArrays.BasisLayout}, Normalized{Float64, LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Normalized{Float64, Legendre{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
    @ ArrayLayouts ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:22
 [13] ldiv
    @ ~\.julia\packages\ArrayLayouts\pPGEV\src\ldiv.jl:86 [inlined]
 [14] \(A::Normalized{Float64, LanczosPolynomial{Float64, QuasiArrays.ApplyQuasiVector{Float64, typeof(*), Tuple{Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, LazyArrays.CachedArray{Float64, 1, Vector{Float64}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}}}, Normalized{Float64, QuasiArrays.SubQuasiArray{Float64, 2, Legendre{Float64}, Tuple{ContinuumArrays.AffineMap{Float64, Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, Inclusion{Float64, ChebyshevInterval{Float64}}}, Base.Slice{InfiniteArrays.OneToInf{Int64}}}, false}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, B::Normalized{Float64, Legendre{Float64}, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}})
    @ QuasiArrays ~\.julia\packages\QuasiArrays\z4zwW\src\matmul.jl:34
 [15] top-level scope
    @ In[30]:6
 [16] eval
    @ .\boot.jl:360 [inlined]
 [17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094

does not work with `Symbols.jl` due to type casting (it seems)

using Symbolics
using ClassicalOrthogonalPolynomials
P = Legendre()
@variables z
P[z,1:10]

breaks with

MethodError: no method matching Float64(::Symbolics.Num)

Closest candidates are:
(::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:760
(::Type{T})(!Matched::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
convert(::Type{Float64}, ::Symbolics.Num)@number.jl:7
to_quasi_index(::Type{Float64}, ::Symbolics.Num)@indices.jl:103
to_quasi_index(::ClassicalOrthogonalPolynomials.Legendre{Float64}, ::Type, ::Symbolics.Num)@indices.jl:111
to_indices@indices.jl:118[inlined]
to_indices@multidimensional.jl:413[inlined]
view@subquasiarray.jl:67[inlined]
layout_getindex@ArrayLayouts.jl:108[inlined]
getindex(::ClassicalOrthogonalPolynomials.Legendre{Float64}, ::Symbolics.Num, ::UnitRange{Int64})@clenshaw.jl:67
top-level scope@Local: 1[inlined]

Ambiguities in ClassicalOrthogonalPolynomials and its dependencies

I've submitted a PR to QuantumOptics.jl which replaces it's own implementation of Hermite polynomials with this package (qojulia/QuantumOpticsBase.jl#137). However adding ClassicalOrthogonalPolynomials as a dependency caused the Aqua tests for ambiguities to break in QuantumOpticsBase due 150 ambiguities total brought in by ClassicalOrthogonalPolynomials and it's dependencies. This can be checked by running

julia> using Test
julia> using ClassicalOrthogonalPolynomials
julia> Test.detect_ambiguities(Base, Core)

The set of packages which produce ambiguities are:

ArrayLayouts
BandedMatrices
BlockArrays
BlockBandedMatrices
ClassicalOrthogonalPolynomials
ContinuumArrays
DomainSets
InfiniteArrays
Infinities
LazyArrays
QuasiArrays
SemiseparableMatrices

Unfortunately there does not appear to be a way to disable ambiguities on a per package basis. It appears all these packages are all under the JuliaApproximation, JuliaArrays, and JuliaMath organizations with overlapping developers so hopefully there's already some awareness of these ambiguities and maybe some plan to eventually fix them?

Laguerre

Is this package expected to, at some point, include the Laguerre polynomials (for evaluation), or are they already in some other package?

`points(Ultraspherical(1, -1..1), n)` returns points on Chebyshev grid of the first kind and not the second kind

using ApproxFunOrthogonalPolynomials

# For reference
chebyshevpoints(Float64, 3, Val(1))  # [0.866, 0.0, -0.866]
chebyshevpoints(Float64, 3, Val(2))  # [1.0, 0.0, -1.0]

# These should be points corresponding to the Chebyshev polynomials of the second kind, I thought?
points(Ultraspherical(1, -1..1), 3)  # [0.866, 0.0, -0.866]

UndefVarError in docs

What's going on in the docs that almost all lines of code evaluate to UndefVarError?

Add Associated(::OrthogonalPolynomial)

I think we want to support e.g. Associated(Legendre()) for the associated OPs. It should be straight forward;

struct Associated{T, PP} <: OrthogonalPolynomial{T}
   P::PP
end

function jacobimatrix(Q::Associated)
    X = jacobimatrix(Q.P)
    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
    Tridiagonal(c[2:end], a[2:end], b[2:end])
end

@MikaelSlevinsky Are there any other features I should add? Do you have that the transforms implemented yet?

We could also support Associated(Legendre(), m) for more general case.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

What's the structure of "OffHilbert"?

@MikaelSlevinsky Any thoughts on how we o]should treat the structure of "OffHilbert"?

At the moment it's just treated as dense, which is not so nice for adaptive linear algebra. I could introduce a notion of "compact".

Quadratically shifted bases

I want to use $P(2x^2-1)$. I had this working when the repo was still OrthogonalPolynomialsQuasi.jl but some changes since have stopped it from working with my setup and I don't want to use the old repo. I think this should just be easy to do out of the box in this package, especially given how ubiquitous $P(2x^2-1)$ is for rotationally symmetric bases.

It seems a Chebyhev quadratic shift is implemented in the tests, see

ClassicalOrthogonalPolynomials.jl/test/test_chebyshev.jl

Line 478 in 139fee6

struct QuadraticMap{T} <: Map{T} end

I think this should be low effort to get working for at least general Jacobi polynomials. Any objections to me adding this functionality to the package itself rather than just have it defined in tests?

Once we do this I would also be happy to type up some documentation on shifted bases. 😁

How to cite the package?

@dlfivefifty What's the preferred way to cite this package as well as others in its ecosystem such as https://github.com/JuliaApproximation/ContinuumArrays.jl and https://github.com/JuliaApproximation/QuasiArrays.jl ? I want to specifically refer to these packages in my thesis, hence the question.

Some considerations:

Any papers associated with them to cite (either along with the repo or instead of)?
Perhaps it would be good to consider either setting up a Zenodo DOI or leaving some notes in the readme file on which papers to cite instead?

copy(Ldiv(A, B)) for Normalized(A) and Weighted(Normalized(B)) ambiguity

using ClassicalOrthogonalPolynomials
P⁰⁰ = Jacobi(0, 0) |> Normalized
P¹⁰ = Jacobi(1, 0) |> Normalized
L₁₀⁰⁰ = P⁰⁰ \ Weighted(P¹⁰)

LoadError: MethodError: copy(::Ldiv{ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}, ClassicalOrthogonalPolynomials.WeightedOPLayout{ClassicalOrthogonalPolynomials.NormalizedOPLayout{ClassicalOrthogonalPolynomials.OPLayout}}, Normalized{Float64, Jacobi{Float64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, Weighted{Float64, Normalized{Float64, Jacobi{Float64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}}) is ambiguous.

Candidates:
  copy(L::Ldiv{<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout, Lay}) where Lay
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\nNR0P\src\normalized.jl:124
  copy(L::Ldiv{<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout, Lay}) where Lay<:ContinuumArrays.AbstractBasisLayout
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\nNR0P\src\normalized.jl:126
  copy(L::Ldiv{Lay, <:ClassicalOrthogonalPolynomials.WeightedOPLayout{<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout}}) where Lay<:ContinuumArrays.AbstractBasisLayout
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\nNR0P\src\normalized.jl:194
  copy(L::Ldiv{<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout, Lay}) where Lay<:LazyArrays.AbstractLazyLayout
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\nNR0P\src\normalized.jl:128
  copy(P::Ldiv{<:ContinuumArrays.AbstractBasisLayout, <:ContinuumArrays.AbstractBasisLayout})
    @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\d6Ead\src\bases\bases.jl:81
  copy(L::Ldiv{<:ContinuumArrays.AbstractBasisLayout, <:LazyArrays.AbstractLazyLayout})
    @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\d6Ead\src\bases\bases.jl:320
  copy(L::Ldiv{<:LazyArrays.AbstractLazyLayout, <:LazyArrays.AbstractLazyLayout})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\nyS8r\src\linalg\inv.jl:124
  copy(L::Ldiv{<:ContinuumArrays.AbstractBasisLayout})
    @ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\d6Ead\src\bases\bases.jl:319
  copy(L::Ldiv{<:Any, <:LazyArrays.AbstractLazyLayout})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\nyS8r\src\linalg\inv.jl:126
  copy(L::Ldiv{<:LazyArrays.AbstractLazyLayout})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\nyS8r\src\linalg\inv.jl:125

Possible fix, define
  copy(::Ldiv{Lay, Lay}) where {Lay<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout, Lay<:(ClassicalOrthogonalPolynomials.WeightedOPLayout{<:ClassicalOrthogonalPolynomials.AbstractNormalizedOPLayout})}

Stacktrace:
  [1] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ldiv.jl:22 [inlined]
  [2] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\feBnP\src\ldiv.jl:98 [inlined]
  [3] \(A::Normalized{Float64, Jacobi{Float64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}, B::Weighted{Float64, Normalized{Float64, Jacobi{Float64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}})
    @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\Jo8rC\src\matmul.jl:34
  [4] top-level scope
    @ c:\Users\User\Documents\PhD\DanielPhDWork\notes\PolynomialsFromFactorisations\infrevchol.jl:127
  [5] eval
    @ .\boot.jl:385 [inlined]
  [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:2076
  [7] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base .\essentials.jl:892
  [8] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:889
  [9] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:271
 [10] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:181
 [11] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [12] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [13] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [14] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:150
 [15] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:515
 [16] with_logger
    @ .\logging.jl:627 [inlined]
 [17] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [18] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [19] invokelatest(::Any)
    @ Base .\essentials.jl:889
 [20] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:34
in expression starting at c:\Users\User\Documents\PhD\DanielPhDWork\notes\PolynomialsFromFactorisations\infrevchol.jl:127

Will make a PR fixing this soon.

Lanczos Jacobi matrices and `^`

Jacobi matrices obtained from Lanczos don't appear to play nice with ^:

julia>  using ClassicalOrthogonalPolynomials

julia> P = Legendre()
Legendre()

julia> x = axes(P,1)
Inclusion(-1.0..1.0 (Chebyshev))

julia> T = LanczosPolynomial(@.(2-x),P)
LanczosPolynomial{Float64} with weight with singularities LegendreWeight{Float64}

julia> X = jacobimatrix(T)
ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos:
 -0.166667  0.552771    ⋅           ⋅           ⋅            ⋅         ⋅    ⋅    ⋅   …  
  0.552771  0.0212121  0.519377     ⋅           ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅        0.519377   0.00679908  0.508021     ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅         0.508021    0.00221202  0.504252      ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅          ⋅          0.504252    0.000936283  0.502645   ⋅    ⋅    ⋅      
  ⋮                                                         ⋮                        ⋱  

julia> X^1
ERROR: MethodError: Cannot `convert` an object of type LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}} to an object of type ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}
Closest candidates are:
  convert(::Type{T}, ::Factorization) where T<:AbstractArray at ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:58
  convert(::Type{T}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1}) where T at ~/.julia/packages/QuasiArrays/tF3vb/src/multidimensional.jl:131
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:16
  ...
Stacktrace:
 [1] to_power_type(x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}})
   @ Base ./intfuncs.jl:250
 [2] power_by_squaring(x_::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
   @ Base ./intfuncs.jl:265
 [3] _power_by_squaring
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:337 [inlined]
 [4] power_by_squaring
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:345 [inlined]
 [5] ^(A::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
   @ LinearAlgebra ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:454
 [6] literal_pow(f::typeof(^), x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, #unused#::Val{1})
   @ Base ./intfuncs.jl:340
 [7] top-level scope
   @ ~/Documents/Projects/ClassicalOrthogonalPolynomials.jl/test/exp.jl:267

julia> X^2
ERROR: MethodError: Cannot `convert` an object of type LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}} to an object of type ApplyArray{Float64, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}
Closest candidates are:
  convert(::Type{T}, ::Factorization) where T<:AbstractArray at ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:58
  convert(::Type{T}, ::QuasiArrays.QuasiIteratorsMD.QuasiCartesianIndex{1}) where T at ~/.julia/packages/QuasiArrays/tF3vb/src/multidimensional.jl:131
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:16
  ...
Stacktrace:
 [1] to_power_type(x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}})
   @ Base ./intfuncs.jl:250
 [2] power_by_squaring(x_::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
   @ Base ./intfuncs.jl:265
 [3] _power_by_squaring
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:337 [inlined]
 [4] power_by_squaring
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:345 [inlined]
 [5] ^(A::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, p::Int64)
   @ LinearAlgebra ~/Documents/Backends/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:454
 [6] literal_pow(f::typeof(^), x::LazyBandedMatrices.SymTridiagonal{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}, #unused#::Val{2})
   @ Base ./intfuncs.jl:340
 [7] top-level scope
   @ ~/Documents/Projects/ClassicalOrthogonalPolynomials.jl/test/exp.jl:268

Multiplication and addition work fine:

julia> X*X
(ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos) * (ℵ₀×ℵ₀ SymTridiagonal{Float64} Jacobi operator from Lanczos) with indices OneToInf()×OneToInf():
  0.333333  -0.080403   0.287096     ⋅           ⋅           ⋅            ⋅         ⋅   …  
 -0.080403   0.575758   0.0145484   0.263854     ⋅           ⋅            ⋅         ⋅      
  0.287096   0.0145484  0.527884    0.00457783  0.256171     ⋅            ⋅         ⋅      
   ⋅         0.263854   0.00457783  0.512361    0.00158754  0.25346       ⋅         ⋅      
   ⋅          ⋅         0.256171    0.00158754  0.506923    0.000714292  0.252233   ⋅      
  ⋮                                                         ⋮                           ⋱  

julia> X+X
ℵ₀×ℵ₀ LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(+), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, BroadcastVector{Float64, typeof(+), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}} with indices OneToInf()×OneToInf():
 -0.333333  1.10554     ⋅          ⋅           ⋅           ⋅        ⋅    ⋅    ⋅    ⋅   …  
  1.10554   0.0424242  1.03875     ⋅           ⋅           ⋅        ⋅    ⋅    ⋅    ⋅      
   ⋅        1.03875    0.0135982  1.01604      ⋅           ⋅        ⋅    ⋅    ⋅    ⋅      
   ⋅         ⋅         1.01604    0.00442404  1.0085       ⋅        ⋅    ⋅    ⋅    ⋅      
   ⋅         ⋅          ⋅         1.0085      0.00187257  1.00529   ⋅    ⋅    ⋅    ⋅      
  ⋮                                                       ⋮                            ⋱

In the meantime, there are many workarounds. E.g.:

julia> X = Symmetric(BandedMatrix(0 => X.dv, 1 => X.ev))
ℵ₀×ℵ₀ Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
 -0.166667  0.552771    ⋅           ⋅           ⋅            ⋅         ⋅    ⋅    ⋅   …  
  0.552771  0.0212121  0.519377     ⋅           ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅        0.519377   0.00679908  0.508021     ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅         0.508021    0.00221202  0.504252      ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅          ⋅          0.504252    0.000936283  0.502645   ⋅    ⋅    ⋅      
  ⋮                                                         ⋮                        ⋱  

julia> X^1
ℵ₀×ℵ₀ Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
 -0.166667  0.552771    ⋅           ⋅           ⋅            ⋅         ⋅    ⋅    ⋅   …  
  0.552771  0.0212121  0.519377     ⋅           ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅        0.519377   0.00679908  0.508021     ⋅            ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅         0.508021    0.00221202  0.504252      ⋅         ⋅    ⋅    ⋅      
   ⋅         ⋅          ⋅          0.504252    0.000936283  0.502645   ⋅    ⋅    ⋅      
  ⋮                                                         ⋮                        ⋱  

julia> X^2
ℵ₀×ℵ₀ Symmetric{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}, Symmetric{Float64, BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(*), Tuple{Matrix{Float64}, ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, InfiniteArrays.ReshapedArray{Float64, 2, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{Int64, Infinities.InfiniteCardinal{0}}, Tuple{}}}}}}}}, InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
  0.333333  -0.080403   0.287096     ⋅           ⋅           ⋅            ⋅         ⋅   …  
 -0.080403   0.575758   0.0145484   0.263854     ⋅           ⋅            ⋅         ⋅      
  0.287096   0.0145484  0.527884    0.00457783  0.256171     ⋅            ⋅         ⋅      
   ⋅         0.263854   0.00457783  0.512361    0.00158754  0.25346       ⋅         ⋅      
   ⋅          ⋅         0.256171    0.00158754  0.506923    0.000714292  0.252233   ⋅      
  ⋮                                                         ⋮                           ⋱

Infinite loop converting between ChebyshevU and Ultraspherical

Noticed accidentally that the code

julia> using ClassicalOrthogonalPolynomials

julia> U = ChebyshevU();

julia> U \ diff(U)

never finishes. Ctrl+Cing this gives

julia> U \ diff(U)
ERROR: InterruptException:
Stacktrace:
  [1] cache_filldata!(::LazyArrays.CachedArray{…}, ::UnitRange{…}, ::Base.OneTo{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:222
  [2] resizedata!(::ArrayLayouts.DenseColumnMajor, ::ArrayLayouts.ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:263
  [3] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:218 [inlined]
  [4] setindex!
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\cache.jl:82 [inlined]
  [5] setindex!
    @ .\subarray.jl:331 [inlined]
  [6] _setindex!
    @ .\abstractarray.jl:1426 [inlined]
  [7] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [8] copyto_unaliased!(deststyle::IndexCartesian, dest::SubArray{…}, srcstyle::IndexCartesian, src::SubArray{…})
    @ Base .\abstractarray.jl:1102
  [9] copyto!
    @ .\abstractarray.jl:1068 [inlined]
 [10] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:255 [inlined]
 [11] _copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:259 [inlined]
 [12] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ArrayLayouts.jl:264 [inlined]
 [13] ldiv!(Y::LazyArrays.CachedArray{…}, A::MatrixFactorizations.QR{…}, B::BandedMatrices.BandedMatrix{…})
    @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.3+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\factorization.jl:179
 [14] _ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:84 [inlined]
 [15] copyto!(dest::LazyArrays.CachedArray{…}, M::ArrayLayouts.Ldiv{…}; kwds::@Kwargs{})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113
 [16] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
 [17] ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:104 [inlined]
 [18] _ldiv!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:83 [inlined]
 [19] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:113 [inlined]
 [20] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:21 [inlined]
 [21] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [22] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [23] copy(M::ArrayLayouts.Mul{…})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:90
 [24] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:127 [inlined]
 [25] mul
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:128 [inlined]
 [26] *(A::LazyArrays.ApplyArray{…}, B::BandedMatrices.BandedMatrix{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\mul.jl:213
 [27] _copy_ldiv_mul
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:112 [inlined]
 [28] copy
    @ C:\Users\User\.julia\packages\LazyArrays\JNdsG\src\linalg\inv.jl:116 [inlined]
 [29] basis_ldiv_size
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:327 [inlined]
 [30] copy
    @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:323 [inlined]
 [31] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [32] ldiv
    @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [33] \(A::ChebyshevU{Float64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
    @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
 [34] top-level scope
    @ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.

I'm guessing that this should give a "Not implemented error" since the equivalent does

julia> Ultraspherical(1) \ diff(U)
ERROR: Not implemented
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] \(C2::Ultraspherical{Float64, Int64}, C1::Ultraspherical{Float64, Int64})
   @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\Hdhud\src\classical\ultraspherical.jl:211
 [3] copy
   @ C:\Users\User\.julia\packages\ContinuumArrays\0grIp\src\bases\bases.jl:0 [inlined]
 [4] materialize
   @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:22 [inlined]
 [5] ldiv
   @ C:\Users\User\.julia\packages\ArrayLayouts\ccyJu\src\ldiv.jl:98 [inlined]
 [6] \(A::Ultraspherical{Float64, Int64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\matmul.jl:34
 [7] top-level scope
   @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

Fourier{BigFloat}

Is there a way to work with Fourier and BigFloat? It goes via FFTW which is causing some issues.

using ClassicalOrthogonalPolynomials
F = Fourier{BigFloat}()
F[:, 1:10] \ sin.(axes(F,1))


ERROR: MethodError: no method matching mul!(::Vector{BigFloat}, ::FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}, ::Vector{BigFloat}, ::Bool, ::Bool)

Closest candidates are:
  mul!(::StridedVecOrMat, ::SparseArrays.AbstractSparseMatrixCSC, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:30
  mul!(::StridedVecOrMat, ::LinearAlgebra.Adjoint{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
  mul!(::StridedVecOrMat, ::LinearAlgebra.Transpose{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
  ...

Stacktrace:
  [1] mul!
    @ C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\matmul.jl:276 [inlined]
  [2] mul!(ret::Vector{BigFloat}, F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, bin::Vector{Float64})
    @ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:128
  [3] *(F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, b::Vector{Float64})
    @ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:133
  [4] \(a::ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\plans.jl:26
  [5] \(a::ContinuumArrays.ProjectionFactorization{BigFloat, ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, UnitRange{Int64}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:213
  [6] transform_ldiv_size(#unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64}, A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:260
  [7] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:261
  [8] basis_ldiv_size
    @ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:310 [inlined]
  [9] copy
    @ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:302 [inlined]
 [10] #materialize#13
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
 [11] materialize
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
 [12] #ldiv#20
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
 [13] ldiv
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
 [14] \(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ QuasiArrays C:\Users\john.papad\.julia\packages\QuasiArrays\GYEdf\src\matmul.jl:34
 [15] top-level scope
    @ Untitled-1:6

long compilation time for Jacobi(m,n) \ Jacobi(0,0)

julia> using ClassicalOrthogonalPolynomials

julia> @time Jacobi(2,2) \ Jacobi(0,0);
  0.206673 seconds (326.39 k allocations: 19.069 MiB, 98.56% compilation time)

julia> @time Jacobi(5,5) \ Jacobi(0,0);
  0.435400 seconds (788.60 k allocations: 43.161 MiB, 98.64% compilation time)

julia> @time Jacobi(10,10) \ Jacobi(0,0);
 16.587771 seconds (26.72 M allocations: 1.142 GiB, 1.86% gc time, 99.94% compilation time)

julia> @time Jacobi(11,11) \ Jacobi(0,0);
 46.804348 seconds (75.86 M allocations: 3.206 GiB, 1.57% gc time, 100.00% compilation time)

julia> @time Jacobi(11,11) \ Jacobi(0,0);
  0.000060 seconds (43 allocations: 20.266 KiB)

This could possibly be caused by LazyArrays.ApplyArray using NTuple.