Support Legendre and Hermite

Author: dlfivefifty

ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl

@@ -1,7 +1,8 @@
 module ClassicalOrthogonalPolynomialsMutableArithmeticsExt
 using ClassicalOrthogonalPolynomials, MutableArithmetics
-using ClassicalOrthogonalPolynomials: initiateforwardrecurrence, recurrencecoefficients, _p0
+import ClassicalOrthogonalPolynomials: initiateforwardrecurrence, recurrencecoefficients, _p0
 
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractMutable, n::Number) = initiateforwardrecurrence(n-1, recurrencecoefficients(P)..., x, _p0(P))[end]
 
+recurrencecoefficients(::Legendre{T}) where T<:AbstractMutable = recurrencecoefficients(Legendre())
 end

src/classical/hermite.jl

@@ -46,15 +46,15 @@ axes(::Hermite{T}) where T = (Inclusion{T}(ℝ), oneto(∞))
 computes the `n`-th Hermite polynomial, orthogonal with 
 respec to `exp(-x^2)`, at `z`.
 """
-hermiteh(n::Integer, z::Number) = Base.unsafe_getindex(Hermite{typeof(z)}(), z, n+1)
+hermiteh(n::Integer, z) = Base.unsafe_getindex(Hermite{typeof(z)}(), z, n+1)
 
 broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), ::HermiteWeight{T}, ::Hermite{V}) where {T,V} = Weighted(Hermite{promote_type(T,V)}())
 
 # H_{n+1} = 2x H_n - 2n H_{n-1}
 # 1/2 * H_{n+1} + n H_{n-1} = x H_n 
 # x*[H_0 H_1 H_2 …] = [H_0 H_1 H_2 …] * [0    1; 1/2  0     2; 1/2   0  3; …]   
 jacobimatrix(H::Hermite{T}) where T = Tridiagonal(Fill(one(T)/2,∞), Zeros{T}(∞), one(T):∞)
-recurrencecoefficients(H::Hermite{T}) where T = Fill{T}(2,∞), Zeros{T}(∞), zero(T):2:∞
+recurrencecoefficients(H::Hermite) where T = Fill(2,∞), Zeros{Int}(∞), 0:2:∞
 
 weightedgrammatrix(::Hermite{T}) where T = Diagonal(sqrt(convert(T,π)) .* convert(T,2) .^ (0:∞) .* gamma.(one(T):∞))
 

src/classical/jacobi.jl

@@ -119,8 +119,8 @@ basis_singularities(w::JacobiWeight) = Weighted(Jacobi(w.a, w.b))
 computes the `n`-th Jacobi polynomial, orthogonal with
 respec to `(1-x)^a*(1+x)^b`, at `z`.
 """
-jacobip(n::Integer, a, b, z::Number) = Base.unsafe_getindex(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
-normalizedjacobip(n::Integer, a, b, z::Number) = Base.unsafe_getindex(Normalized(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b)), z, n+1)
+jacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
+normalizedjacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Normalized(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b)), z, n+1)
 
 OrthogonalPolynomial(w::JacobiWeight) = Jacobi(w.a, w.b)
 orthogonalityweight(P::Jacobi) = JacobiWeight(P.a, P.b)

src/classical/ultraspherical.jl

@@ -54,7 +54,7 @@ const WeightedUltraspherical{T} = WeightedBasis{T,<:UltrasphericalWeight,<:Ultra
 
 orthogonalityweight(C::Ultraspherical) = UltrasphericalWeight(C.λ)
 
-ultrasphericalc(n::Integer, λ, z::Number) = Base.unsafe_getindex(Ultraspherical{promote_type(typeof(λ),typeof(z))}(λ), z, n+1)
+ultrasphericalc(n::Integer, λ, z) = Base.unsafe_getindex(Ultraspherical{promote_type(typeof(λ),typeof(z))}(λ), z, n+1)
 ultraspherical(λ, d::AbstractInterval{T}) where T = Ultraspherical{float(promote_type(eltype(λ),T))}(λ)[affine(d,ChebyshevInterval{T}()), :]
 
 ==(a::Ultraspherical, b::Ultraspherical) = a.λ == b.λ

test/test_dynamicpolynomials.jl

@@ -6,6 +6,7 @@ using DynamicPolynomials, ClassicalOrthogonalPolynomials, Test
     @test chebyshevt(0,x) == 1
     @test chebyshevt(1,x) == x
     @test chebyshevt(5,x) == 5x - 20x^3 + 16x^5
-    @test chebyshevu(5,x) == 6x - 32x^3 + 32x^5
-    legendrep(5,x)
+    @test chebyshevu(5,x) == ultrasphericalc(5,1,x) == 6x - 32x^3 + 32x^5
+    @test legendrep(5,x) ≈ (15x - 70x^3 + 63x^5)/8
+    @test hermiteh(5,x) == 120x - 160x^3 + 32x^5
 end