Support DynamicPolynomials.jl (#214)

* Support DynamicPolynomials.jl

* Support Legendre and Hermite

* Update hermite.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.13.7"
+version = "0.13.8"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -26,13 +26,20 @@ RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 RecurrenceRelationshipArrays = "b889d2dc-af3c-4820-88a8-238fa91d3518"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
+[weakdeps]
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+
+[extensions]
+ClassicalOrthogonalPolynomialsMutableArithmeticsExt = "MutableArithmetics"
+
 [compat]
 ArrayLayouts = "1.3.1"
 BandedMatrices = "1"
 BlockArrays = "1"
 BlockBandedMatrices = "0.13"
 ContinuumArrays = "0.18.3"
 DomainSets = "0.6, 0.7"
+DynamicPolynomials = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "1"
 FastTransforms = "0.16.6"
@@ -43,6 +50,7 @@ InfiniteLinearAlgebra = "0.8"
 IntervalSets = "0.7"
 LazyArrays = "2.2"
 LazyBandedMatrices = "0.10"
+MutableArithmetics = "1"
 QuasiArrays = "0.11"
 RecurrenceRelationships = "0.1"
 RecurrenceRelationshipArrays = "0.1"
@@ -51,11 +59,12 @@ julia = "1.10"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random"]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random", "DynamicPolynomials"]

ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl

@@ -0,0 +1,8 @@
+module ClassicalOrthogonalPolynomialsMutableArithmeticsExt
+using ClassicalOrthogonalPolynomials, MutableArithmetics
+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/chebyshev.jl

@@ -61,13 +61,13 @@ chebyshevu(S::AbstractQuasiMatrix) = chebyshevu(axes(S,1))
 
 computes the `n`-th Chebyshev polynomial of the first kind at `z`.
 """
-chebyshevt(n::Integer, z::Number) = Base.unsafe_getindex(ChebyshevT{typeof(z)}(), z, n+1)
+chebyshevt(n::Integer, z) = Base.unsafe_getindex(ChebyshevT{typeof(z)}(), z, n+1)
 """
      chebyshevt(n, z)
 
 computes the `n`-th Chebyshev polynomial of the second kind at `z`.
 """
-chebyshevu(n::Integer, z::Number) = Base.unsafe_getindex(ChebyshevU{typeof(z)}(), z, n+1)
+chebyshevu(n::Integer, z) = Base.unsafe_getindex(ChebyshevU{typeof(z)}(), z, n+1)
 
 chebysevtweight(d::AbstractInterval{T}) where T = ChebyshevTWeight{float(T)}[affine(d,ChebyshevInterval{T}())]
 chebysevuweight(d::AbstractInterval{T}) where T = ChebyshevUWeight{float(T)}[affine(d,ChebyshevInterval{T}())]

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) = 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/legendre.jl

@@ -67,7 +67,7 @@ legendre(d::Inclusion) = legendre(d.domain)
 
 computes the `n`-th Legendre polynomial at `z`.
 """
-legendrep(n::Integer, z::Number) = Base.unsafe_getindex(Legendre{typeof(z)}(), z, n+1)
+legendrep(n::Integer, z) = Base.unsafe_getindex(Legendre{typeof(z)}(), z, n+1)
 
 
 show(io::IO, w::Legendre{Float64}) = summary(io, w)

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/runtests.jl

@@ -85,4 +85,6 @@ end
 @testset "Issue #179" begin
     @test startswith(sprint(show, MIME"text/plain"(), Chebyshev()[0.3, :]; context=(:compact=>true, :limit=>true)), "ℵ₀-element view(::ChebyshevT{Float64}, 0.3, :)")
     @test startswith(sprint(show, MIME"text/plain"(), Jacobi(0.2, 0.5)[-0.7, :]; context=(:compact=>true, :limit=>true)), "ℵ₀-element view(::Jacobi{Float64}, -0.7, :)")
-end
+end
+
+include("test_dynamicpolynomials.jl")

test/test_dynamicpolynomials.jl

@@ -0,0 +1,12 @@
+using DynamicPolynomials, ClassicalOrthogonalPolynomials, Test
+
+@polyvar x
+
+@testset "DynamicPolynomials" begin
+    @test chebyshevt(0,x) == 1
+    @test chebyshevt(1,x) == x
+    @test chebyshevt(5,x) == 5x - 20x^3 + 16x^5
+    @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