Separate type of parameters in Jacobi (#220)

* Update Project.toml

* Update test_dynamicpolynomials.jl

* Separate type of parameters in Jacobi

* Import polynomialtype

* fix tests

* Update jacobi.jl

* Update jacobi.jl

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.13.9"
+version = "0.14.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -52,8 +52,8 @@ LazyArrays = "2.2"
 LazyBandedMatrices = "0.10, 0.11"
 MutableArithmetics = "1"
 QuasiArrays = "0.11"
-RecurrenceRelationshipArrays = "0.1"
-RecurrenceRelationships = "0.1"
+RecurrenceRelationshipArrays = "0.1.2"
+RecurrenceRelationships = "0.1.1"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -41,7 +41,7 @@ import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, Abstract
     plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
 import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg, pochhammer
 import RecurrenceRelationships: forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!,
-                        check_clenshaw_recurrences
+                        check_clenshaw_recurrences, polynomialtype
 import RecurrenceRelationshipArrays: initiateforwardrecurrence, Clenshaw
 import FastGaussQuadrature: jacobimoment
 
@@ -254,7 +254,6 @@ function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof
 end
 
 
-
 # function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::BroadcastQuasiVector, C::OrthogonalPolynomial)
 #     axes(a,1) == axes(C,1) || throw(DimensionMismatch())
 #     # re-expand in OP basis

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) = Base.unsafe_getindex(ChebyshevT{typeof(z)}(), z, n+1)
+chebyshevt(n::Integer, z) = Base.unsafe_getindex(ChebyshevT{polynomialtype(typeof(z))}(), z, n+1)
 """
      chebyshevt(n, z)
 
 computes the `n`-th Chebyshev polynomial of the second kind at `z`.
 """
-chebyshevu(n::Integer, z) = Base.unsafe_getindex(ChebyshevU{typeof(z)}(), z, n+1)
+chebyshevu(n::Integer, z) = Base.unsafe_getindex(ChebyshevU{polynomialtype(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,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) = Base.unsafe_getindex(Hermite{typeof(z)}(), z, n+1)
+hermiteh(n::Integer, z) = Base.unsafe_getindex(Hermite{polynomialtype(typeof(z))}(), z, n+1)
 
 broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), ::HermiteWeight{T}, ::Hermite{V}) where {T,V} = Weighted(Hermite{promote_type(T,V)}())
 

src/classical/jacobi.jl

@@ -136,14 +136,14 @@ The quasi-matrix representing Jacobi polynomials, where the first axes represent
 The eltype, when not specified, will be converted to a floating point data type.
 # Examples
 ```jldoctest
-julia> J=Jacobi(0,0) # The eltype will be converted to float
-Jacobi(0.0, 0.0)
+julia> J=Jacobi(0, 0) # The eltype will be converted to float
+Jacobi(0, 0)
 
 julia> axes(J)
 (Inclusion(-1.0 .. 1.0 (Chebyshev)), OneToInf())
 
 julia> J[0,:] # Values of polynomials at x=0
-ℵ₀-element view(::Jacobi{Float64}, 0.0, :) with eltype Float64 with indices OneToInf():
+ℵ₀-element view(::Jacobi{Float64, $Int}, 0.0, :) with eltype Float64 with indices OneToInf():
   1.0
   0.0
  -0.5
@@ -162,12 +162,14 @@ julia> J0[0],J0[0.5]
 (1.0, 1.0)
 ```
 """
-struct Jacobi{T} <: AbstractJacobi{T}
-    a::T
-    b::T
-    Jacobi{T}(a, b) where T = new{T}(convert(T,a), convert(T,b))
+struct Jacobi{T,V} <: AbstractJacobi{T}
+    a::V
+    b::V
+    Jacobi{T,V}(a, b) where {T,V} = new{T,V}(convert(V,a), convert(V,b))
 end
 
+Jacobi{T}(a::V, b::V) where {T,V} = Jacobi{T,V}(a, b)
+Jacobi{T}(a, b) where T = Jacobi{T}(promote(a,b)...)
 Jacobi(a::V, b::T) where {T,V} = Jacobi{float(promote_type(T,V))}(a, b)
 
 AbstractQuasiArray{T}(w::Jacobi) where T = Jacobi{T}(w.a, w.b)
@@ -181,13 +183,13 @@ The [`Jacobi`](@ref) polynomials affine-mapped to interval `d`.
 # Examples
 ```jldoctest
 julia> J = jacobi(1, 1, 0..1)
-Jacobi(1.0, 1.0) affine mapped to 0 .. 1
+Jacobi(1, 1) affine mapped to 0 .. 1
 
 julia> axes(J)
 (Inclusion(0 .. 1), OneToInf())
 
 julia> J[0,:]
-ℵ₀-element view(::Jacobi{Float64}, -1.0, :) with eltype Float64 with indices OneToInf():
+ℵ₀-element view(::Jacobi{Float64, $Int}, -1.0, :) with eltype Float64 with indices OneToInf():
    1.0
   -2.0
    3.0
@@ -215,8 +217,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) = 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)
+jacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Jacobi{polynomialtype(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
+normalizedjacobip(n::Integer, a, b, z) = Base.unsafe_getindex(Normalized(Jacobi{polynomialtype(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)
@@ -273,15 +275,16 @@ axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()),
 ==(P::Legendre, Q::Jacobi) = Jacobi(P) == Q
 ==(P::Jacobi, Q::Legendre) = P == Jacobi(Q)
 ==(A::WeightedJacobi, B::WeightedJacobi) = A.args == B.args
-==(A::WeightedJacobi, B::Jacobi{T}) where T = A == JacobiWeight(zero(T),zero(T)).*B
+==(A::WeightedJacobi, B::Jacobi{T,V}) where {T,V} = A == JacobiWeight(zero(V),zero(V)).*B
 ==(A::WeightedJacobi, B::Legendre) = A == Jacobi(B)
 ==(A::Legendre, B::WeightedJacobi) = Jacobi(A) == B
-==(A::Jacobi{T}, B::WeightedJacobi) where T = JacobiWeight(zero(T),zero(T)).*A == B
+==(A::Jacobi{T,V}, B::WeightedJacobi) where {T,V} = JacobiWeight(zero(V),zero(V)).*A == B
 ==(A::Legendre, B::Weighted{<:Any,<:AbstractJacobi}) = A == B.P
 ==(A::Weighted{<:Any,<:AbstractJacobi}, B::Legendre) = A.P == B
 
 show(io::IO, P::Jacobi) = summary(io, P)
-summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
+summary(io::IO, P::Jacobi{Float64}) = print(io, "Jacobi($(P.a), $(P.b))")
+summary(io::IO, P::Jacobi{T}) where T = print(io, "Jacobi{$T}($(P.a), $(P.b))")
 
 ###
 # transforms
@@ -523,22 +526,22 @@ diff(S::Jacobi; dims=1) = ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), _BandedMatrix
 
 function diff(S::Jacobi{T}, m::Integer; dims=1) where T
     D = _BandedMatrix((pochhammer.((S.a + S.b+1):∞, m)/convert(T, 2)^m)', ℵ₀, -m, m)
-    ApplyQuasiMatrix(*, Jacobi(S.a+m,S.b+m), D)
+    ApplyQuasiMatrix(*, Jacobi{T}(S.a+m,S.b+m), D)
 end
 
 
 #L_6^t
-function diff(WS::HalfWeighted{:a,<:Any,<:Jacobi}; dims=1)
+function diff(WS::HalfWeighted{:a,T,<:Jacobi}; dims=1) where T
     S = WS.P
     a,b = S.a, S.b
-    ApplyQuasiMatrix(*, HalfWeighted{:a}(Jacobi(a-1,b+1)), Diagonal(-(a:∞)))
+    ApplyQuasiMatrix(*, HalfWeighted{:a}(Jacobi{T}(a-1,b+1)), Diagonal(-(a:∞)))
 end
 
 #L_6
-function diff(WS::HalfWeighted{:b,<:Any,<:Jacobi}; dims=1)
+function diff(WS::HalfWeighted{:b,T,<:Jacobi}; dims=1) where T
     S = WS.P
     a,b = S.a, S.b
-    ApplyQuasiMatrix(*, HalfWeighted{:b}(Jacobi(a+1,b-1)), Diagonal(b:∞))
+    ApplyQuasiMatrix(*, HalfWeighted{:b}(Jacobi{T}(a+1,b-1)), Diagonal(b:∞))
 end
 
 for ab in (:(:a), :(:b))
@@ -549,7 +552,7 @@ for ab in (:(:a), :(:b))
 end
 
 
-function diff(WS::Weighted{<:Any,<:Jacobi}; dims=1)
+function diff(WS::Weighted{T,<:Jacobi}; dims=1) where T
     # L_1^t
     S = WS.P
     a,b = S.a, S.b
@@ -560,13 +563,13 @@ function diff(WS::Weighted{<:Any,<:Jacobi}; dims=1)
     elseif iszero(b)
         diff(HalfWeighted{:a}(S))
     else
-        ApplyQuasiMatrix(*, Weighted(Jacobi(a-1, b-1)), _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1))
+        ApplyQuasiMatrix(*, Weighted(Jacobi{T}(a-1, b-1)), _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1))
     end
 end
 
 
 # Jacobi(a-1,b-1)\ (D*w*Jacobi(a,b))
-function diff(WS::WeightedJacobi; dims=1)
+function diff(WS::WeightedJacobi{T}; dims=1) where T
     w,S = WS.args
     a,b = S.a, S.b
     if isorthogonalityweighted(WS) # L_1^t
@@ -582,22 +585,22 @@ function diff(WS::WeightedJacobi; dims=1)
         # D * ((1+x)^w.b * P^(a,b)) == D * ((1+x)^(w.b-b) * (1+x)^b * P^(a,b))
         #    == (1+x)^(w.b-1) * (w.b-b) * P^(a,b) + (1+x)^(w.b-b) * D*((1+x)^b*P^(a,b))
         #    == (1+x)^(w.b-1) * P^(a+1,b) ((w.b-b) * C2 + C1 * W)
-        W = HalfWeighted{:b}(Jacobi(a+1, b-1)) \ diff(HalfWeighted{:b}(S))
-        J = Jacobi(a+1,b) # range Jacobi
-        C1 = J \ Jacobi(a+1, b-1)
-        C2 = J \ Jacobi(a,b)
+        W = HalfWeighted{:b}(Jacobi{T}(a+1, b-1)) \ diff(HalfWeighted{:b}(S))
+        J = Jacobi{T}(a+1,b) # range Jacobi
+        C1 = J \ Jacobi{T}(a+1, b-1)
+        C2 = J \ Jacobi{T}(a,b)
         ApplyQuasiMatrix(*, JacobiWeight(w.a,w.b-1) .* J, (w.b-b) * C2 + C1 * W)
     elseif iszero(w.b)
-        W = HalfWeighted{:a}(Jacobi(a-1, b+1)) \ diff(HalfWeighted{:a}(S))
-        J = Jacobi(a,b+1) # range Jacobi
-        C1 = J \ Jacobi(a-1, b+1)
-        C2 = J \ Jacobi(a,b)
+        W = HalfWeighted{:a}(Jacobi{T}(a-1, b+1)) \ diff(HalfWeighted{:a}(S))
+        J = Jacobi{T}(a,b+1) # range Jacobi
+        C1 = J \ Jacobi{T}(a-1, b+1)
+        C2 = J \ Jacobi{T}(a,b)
         ApplyQuasiMatrix(*, JacobiWeight(w.a-1,w.b) .* J, -(w.a-a) * C2 + C1 * W)
     elseif iszero(a) && iszero(b) # Legendre
         # D * ((1+x)^w.b * (1-x)^w.a * P))
         #    == (1+x)^(w.b-1) * (1-x)^(w.a-1) * ((1-x) * (w.b) * P - (1+x) * w.a * P + (1-x^2) * D * P)
         #    == (1+x)^(w.b-1) * (1-x)^(w.a-1) * ((1-x) * (w.b) * P - (1+x) * w.a * P + P * L * W)
-        J = Jacobi(a+1,b+1) # range space
+        J = Jacobi{T}(a+1,b+1) # range space
         W = J \ diff(S)
         X = jacobimatrix(S)
         L = S \ Weighted(J)
@@ -608,9 +611,9 @@ function diff(WS::WeightedJacobi; dims=1)
         #    == (1+x)^(w.b-1) * (1-x)^(w.a-1) * ((1-x) * (w.b-b) * P^(a,b) + (1+x) * (a-w.a) * P^(a,b))
         #        + (1+x)^(w.b-b) * (1-x)^(w.a-a) * D * ((1+x)^b * (1-x)^a * P^(a,b)))
 
-        W = Weighted(Jacobi(a-1,b-1)) \ diff(Weighted(S))
+        W = Weighted(Jacobi{T}(a-1,b-1)) \ diff(Weighted(S))
         X = jacobimatrix(S)
-        C = S \ Jacobi(a-1,b-1)
+        C = S \ Jacobi{T}(a-1,b-1)
         (JacobiWeight(w.a-1,w.b-1) .* S) *  (((w.b-b+a-w.a)*I+(a-w.a-w.b+b) * X) + C*W)
     end
 end

src/classical/laguerre.jl

@@ -45,7 +45,7 @@ axes(::Laguerre{T}) where T = (Inclusion(HalfLine{T}()), oneto(∞))
 computes the `n`-th generalized Laguerre polynomial, orthogonal with 
 respec to `x^α * exp(-x)`, at `z`.
 """
-laguerrel(n::Integer, α, z::Number) = Base.unsafe_getindex(Laguerre{promote_type(typeof(α), typeof(z))}(α), z, n+1)
+laguerrel(n::Integer, α, z::Number) = Base.unsafe_getindex(Laguerre{polynomialtype(typeof(α), typeof(z))}(α), z, n+1)
 
 """
      laguerrel(n, z)

src/classical/ultraspherical.jl

@@ -54,7 +54,7 @@ const WeightedUltraspherical{T} = WeightedBasis{T,<:UltrasphericalWeight,<:Ultra
 
 orthogonalityweight(C::Ultraspherical) = UltrasphericalWeight(C.λ)
 
-ultrasphericalc(n::Integer, λ, z) = Base.unsafe_getindex(Ultraspherical{promote_type(typeof(λ),typeof(z))}(λ), z, n+1)
+ultrasphericalc(n::Integer, λ, z) = Base.unsafe_getindex(Ultraspherical{polynomialtype(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

@@ -84,7 +84,7 @@ 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, :)")
+    @test startswith(sprint(show, MIME"text/plain"(), Jacobi(0.2, 0.5)[-0.7, :]; context=(:compact=>true, :limit=>true)), "ℵ₀-element view(::Jacobi{Float64, Float64}, -0.7, :)")
 end
 
 include("test_dynamicpolynomials.jl")

test/test_dynamicpolynomials.jl

@@ -9,4 +9,6 @@ using DynamicPolynomials, ClassicalOrthogonalPolynomials, Test
     @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
+    @test jacobip(5,1,0,x) ≈ (5 + 35x - 70x^2 - 210x^3 + 105x^4 + 231x^5)/16
+    @test jacobip(5,0.1,-0.2,x) ≈ 3*(3306827 + 75392855x - 37466770x^2 - 350553810x^3 + 47705335x^4 + 314855211x^5)/128000000
 end