diff --git a/Project.toml b/Project.toml index b122c77..11d58d8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.13.7" +version = "0.13.8" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" @@ -26,6 +26,12 @@ 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" @@ -33,6 +39,7 @@ 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,6 +59,7 @@ 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" @@ -58,4 +67,4 @@ 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"] diff --git a/ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl b/ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl new file mode 100644 index 0000000..10c894f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/classical/chebyshev.jl b/src/classical/chebyshev.jl index 329af2f..210e71e 100644 --- a/src/classical/chebyshev.jl +++ b/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}())] diff --git a/src/classical/hermite.jl b/src/classical/hermite.jl index 0ce139f..911bff8 100644 --- a/src/classical/hermite.jl +++ b/src/classical/hermite.jl @@ -46,7 +46,7 @@ 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)}()) @@ -54,7 +54,7 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), ::HermiteWeight{T}, ::Hermite # 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):∞)) diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl index 2aab7f4..88b3c42 100644 --- a/src/classical/jacobi.jl +++ b/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) diff --git a/src/classical/legendre.jl b/src/classical/legendre.jl index 91a6685..4b4ba2c 100644 --- a/src/classical/legendre.jl +++ 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) diff --git a/src/classical/ultraspherical.jl b/src/classical/ultraspherical.jl index 5cd6df3..6376778 100644 --- a/src/classical/ultraspherical.jl +++ 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.λ diff --git a/test/runtests.jl b/test/runtests.jl index 8148153..4d1190e 100644 --- a/test/runtests.jl +++ 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 \ No newline at end of file +end + +include("test_dynamicpolynomials.jl") \ No newline at end of file diff --git a/test/test_dynamicpolynomials.jl b/test/test_dynamicpolynomials.jl new file mode 100644 index 0000000..5f406b9 --- /dev/null +++ b/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 \ No newline at end of file