using ClassicalOrthogonalPolynomials
using SingularIntegrals
import ClassicalOrthogonalPolynomials: recurrenecoefficients
import SingularIntegrals: RecurrenceArray
P = Jacobi(1/3,1/3)
xc = range(-1,1; length=10000)
n = 10000
# getindex baseline ~ 1.102867 seconds
@time pp = P[xc,1:n];
# Lazy is fast ~ 0.000466 seconds
@time pr = RecurrenceArray(xc, recurrencecoefficients(P), P[xc, 1:2]');
# Actually getting the matrix is "slow" ~ 1.706700 seconds
@time pr[1:n,1:length(xc)]';
pr ≈ pp
Sometimes
RecurrenceArrayis slower thangetindexwhen actually assembling the values e.g.