WIP: Explicit basis

.github/workflows/downstream.yml

@@ -14,12 +14,15 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia-version: [1,1.6]
+        julia-version: [1]
         os: [ubuntu-latest]
         package:
-          - {user: jverzani, repo: SpecialPolynomials.jl, group: All}
           - {user: JuliaControl, repo: ControlSystems.jl, group: All}
-
+          - {user: andreasvarga, repo: DescriptorSystems.jl, group: All}
+          - {user: JuliaDSP,     repo: DSP.jl, group: All}
+          - {user: tkluck,       repo: GaloisFields.jl, group: All}
+          - {user: jverzani,     repo: SpecialPolynomials.jl, group: All}
+          - {user: JuliaGNI,     repo: QuadratureRules.jl, group: All}
     steps:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1

.gitignore

@@ -5,3 +5,4 @@ docs/site/
 *.jl.mem
 
 Manifest.toml
+archive/

Project.toml

@@ -10,6 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -26,16 +27,17 @@ ChainRulesCore = "1"
 MakieCore = "0.6"
 MutableArithmetics = "1"
 RecipesBase = "0.7, 0.8, 1"
+Setfield = "1"
 julia = "1.6"
 
 [extras]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/extending.md

@@ -1,6 +1,6 @@
 # Extending Polynomials
 
-The [`AbstractPolynomial`](@ref) type was made to be extended via a rich interface.
+The [`AbstractPolynomial`](@ref) type was made to be extended via a rich interface; examples follow. The newer [`AbstractUnivariatePolynomial`](@ref) type is illustrated at the end.
 
 ```@docs
 AbstractPolynomial
@@ -148,3 +148,190 @@ julia> p .+ 2
 ```
 
 The unexported `Polynomials.PnPolynomial` type implements much of this.
+
+
+## Extending the AbstractUnivariatePolynomial type
+
+An `AbstractUnivariatePolynomial` polynomial consists of a basis and a storage type. The storage type can be mutable dense, mutable sparse, or immutable dense.
+
+A basis inherits from `Polynomials.AbstractBasis`, in the example our basis type has a parameter. The `ChebyshevT` type, gives a related example of how this task can be implemented.
+
+### The generalized Laguerre polynomials
+
+These are orthogonal polynomials parameterized  by $\alpha$ and defined recursively by
+
+```math
+\begin{align*}
+L^\alpha_1(x) &= 1\\
+L^\alpha_2(x) &= 1 + \alpha - x\\
+L^\alpha_{n+1}(x) &= \frac{2n+1+\alpha -x}{n+1} L^\alpha_n(x) - \frac{n+\alpha}{n+1} L^\alpha_{n-1}(x)\\
+&= (A_nx +B_n) \cdot L^\alpha_n(x) - C_n \cdot L^\alpha_{n-1}(x).
+\end{align*}
+```
+
+There are other [characterizations available](https://en.wikipedia.org/wiki/Laguerre_polynomials). The three-point recursion, described by `A`,`B`, and `C` is used below for evaluation.
+
+We define the basis with
+
+```jldoctest abstract_univariate_polynomial
+julia> using Polynomials;
+
+julia> import Polynomials: AbstractUnivariatePolynomial, AbstractBasis, MutableDensePolynomial;
+
+julia> struct LaguerreBasis{alpha} <: AbstractBasis end
+
+julia> Polynomials.basis_symbol(::Type{<:AbstractUnivariatePolynomial{LaguerreBasis{α},T,X}}) where {α,T,X} =
+           "L^$(α)"
+```
+
+The basis symbol has a poor default. The method requires the full type, as the indeterminate, `X`, may part of the desired output. We added a method to `basis_symbol` to show this basis. More generally, `Polynomials.printbasis` can have methods added to adjust for different display types.
+
+Polynomials can be initiated through specifying a storage type and a basis, say:
+
+```jldoctest abstract_univariate_polynomial
+julia> P = MutableDensePolynomial{LaguerreBasis{0}}
+MutableDensePolynomial{LaguerreBasis{0}}
+
+julia> p = P([1,2,3])
+MutableDensePolynomial(1L^0_0 + 2*L^0_1 + 3*L^0_2)
+```
+
+Or using other storage types:
+
+```jldoctest abstract_univariate_polynomial
+julia> Polynomials.ImmutableDensePolynomial{LaguerreBasis{1}}((1,2,3))
+Polynomials.ImmutableDensePolynomial(1L^1_0 + 2*L^1_1 + 3*L^1_2)
+```
+
+All polynomials have vector addition and scalar multiplication defined:
+
+```jldoctest abstract_univariate_polynomial
+julia> q = P([1,2])
+MutableDensePolynomial(1L^0_0 + 2*L^0_1)
+
+julia> p + q
+MutableDensePolynomial(2L^0_0 + 4*L^0_1 + 3*L^0_2)
+```
+
+```jldoctest abstract_univariate_polynomial
+julia> 2p
+MutableDensePolynomial(2L^0_0 + 4*L^0_1 + 6*L^0_2)
+```
+
+For a new basis, there are no default methods for polynomial evaluation, scalar addition, and polynomial multiplication; and no defaults for `one`, and `variable`.
+
+For the Laguerre Polynomials, Clenshaw recursion can be used for evaluation. Internally, `evalpoly` is called so we forward that method.
+
+```jldoctest abstract_univariate_polynomial
+julia> function ABC(::Type{LaguerreBasis{α}}, n) where {α}
+           o = one(α)
+           d = n + o
+           (A=-o/d, B=(2n + o + α)/d, C=(n+α)/d)
+       end
+ABC (generic function with 1 method)
+```
+
+```jldoctest abstract_univariate_polynomial
+julia> function clenshaw_eval(p::P, x::S) where {α, Bᵅ<: LaguerreBasis{α}, T, P<:AbstractUnivariatePolynomial{Bᵅ,T}, S}
+           d = degree(p)
+           R = typeof(((one(α) * one(T)) * one(S)) / 1)
+           p₀ = one(R)
+           d == -1 && return zero(R)
+           d == 0 && return p[0] * one(R)
+           Δ0 = p[d-1]
+           Δ1 = p[d]
+           @inbounds for i in (d - 1):-1:1
+               A,B,C = ABC(Bᵅ, i)
+               Δ0, Δ1 =
+                   p[i] - Δ1 * C, Δ0 + Δ1 * muladd(x, A, B)
+           end
+           A,B,C = ABC(Bᵅ, 0)
+           p₁ = muladd(x, A, B) * p₀
+           return Δ0 * p₀ + Δ1 * p₁
+       end
+clenshaw_eval (generic function with 1 method)
+```
+
+```jldoctest abstract_univariate_polynomial
+julia> Polynomials.evalpoly(x, p::P) where {P<:AbstractUnivariatePolynomial{<:LaguerreBasis}} =
+               clenshaw_eval(p, x)
+```
+
+We test it out by passing in the variable `x` in the standard basis:
+
+```jldoctest abstract_univariate_polynomial
+julia> p = P([0,0,1])
+MutableDensePolynomial(L^0_2)
+
+julia> x = variable(Polynomial)
+Polynomial(1.0*x)
+
+julia> p(x)
+Polynomial(1.0 - 2.0*x + 0.5*x^2)
+```
+
+We see that conversion to the `Polynomial` type is available through polynomial evaluation. This is used by default, so we have `convert` methods available:
+
+```jldoctest abstract_univariate_polynomial
+julia> convert(ChebyshevT, p)
+ChebyshevT(1.25⋅T_0(x) - 2.0⋅T_1(x) + 0.25⋅T_2(x))
+```
+
+Or, using some extra annotations to have rational arithmetic used, we can compare to easily found representations in the standard basis:
+
+```jldoctest abstract_univariate_polynomial
+julia> q = Polynomials.basis(MutableDensePolynomial{LaguerreBasis{0//1}, Int}, 5)
+MutableDensePolynomial(L^0//1_5)
+
+julia> x = variable(Polynomial{Int})
+Polynomial(x)
+
+julia> q(x)
+Polynomial(1//1 - 5//1*x + 5//1*x^2 - 5//3*x^3 + 5//24*x^4 - 1//120*x^5)
+```
+
+
+To implement scalar addition, we utilize the fact that ``L_0 = 1`` to manipulate the coefficients. Below we specialize to a container type:
+
+```jldoctest abstract_univariate_polynomial
+julia> function Polynomials.scalar_add(c::S, p::P) where {B<:LaguerreBasis,T,X,
+                                                          P<:MutableDensePolynomial{B,T,X},S}
+           R = promote_type(T,S)
+           iszero(p) && return MutableDensePolynomial{B,R,X}(c)
+           cs = convert(Vector{R}, copy(p.coeffs))
+           cs[1] += c
+           MutableDensePolynomial{B,R,X}(cs)
+       end
+
+julia> p + 3
+MutableDensePolynomial(3L^0_0 + L^0_2)
+```
+
+The values of `one` and `variable` are straightforward, as ``L_0=1`` and ``L_1=1 - x`` or ``x = 1 - L_1``
+
+```jldoctest abstract_univariate_polynomial
+julia> Polynomials.one(::Type{P}) where {B<:LaguerreBasis,T,X,P<:AbstractUnivariatePolynomial{B,T,X}} =
+           P([one(T)])
+
+julia> Polynomials.variable(::Type{P}) where {B<:LaguerreBasis,T,X,P<:AbstractUnivariatePolynomial{B,T,X}} =
+           P([one(T), -one(T)])
+```
+
+To see this is correct, we have:
+
+```jldoctest abstract_univariate_polynomial
+julia> variable(P)(x) == x
+true
+```
+
+Finally, we implement polynomial multiplication through conversion to the polynomial type. The [direct formula](https://londmathsoc.onlinelibrary.wiley.com/doi/pdf/10.1112/jlms/s1-36.1.399) could be implemented.
+
+```jldoctest abstract_univariate_polynomial
+julia> function Base.:*(p::MutableDensePolynomial{B,T,X},
+                        q::MutableDensePolynomial{B,S,X}) where {B<:LaguerreBasis, T,S,X}
+           x = variable(Polynomial{T,X})
+           p(x) * q(x)
+       end
+```
+
+Were it defined, a `convert` method from `Polynomial` to the `LaguerreBasis` could be used to implement multiplication.

docs/src/index.md

@@ -790,6 +790,7 @@ julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues`
            end
        end
 
+
 julia> d
 ((x - 4.0) * (x - 1.0000000000000002)) // ((x - 5.0) * (x - 2.0))
 ```

src/Polynomials.jl

@@ -3,6 +3,7 @@ module Polynomials
 #  using GenericLinearAlgebra ## remove for now. cf: https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/71#issuecomment-743928205
 using LinearAlgebra
 import Base: evalpoly
+using Setfield
 
 include("abstract.jl")
 include("show.jl")
@@ -15,14 +16,30 @@ include("common.jl")
 # Polynomials
 include("polynomials/standard-basis.jl")
 include("polynomials/Polynomial.jl")
-include("polynomials/ImmutablePolynomial.jl")
-include("polynomials/SparsePolynomial.jl")
-include("polynomials/LaurentPolynomial.jl")
-include("polynomials/pi_n_polynomial.jl")
 include("polynomials/factored_polynomial.jl")
+
+# polynomials with explicit basis
+include("abstract-polynomial.jl")
+include("polynomial-basetypes/mutable-dense-polynomial.jl")
+include("polynomial-basetypes/mutable-dense-view-polynomial.jl")
+include("polynomial-basetypes/mutable-dense-laurent-polynomial.jl")
+include("polynomial-basetypes/immutable-dense-polynomial.jl")
+include("polynomial-basetypes/mutable-sparse-polynomial.jl")
+
+include("standard-basis/standard-basis.jl")
+include("standard-basis/polynomial.jl")
+include("standard-basis/pn-polynomial.jl")
+include("standard-basis/laurent-polynomial.jl")
+include("standard-basis/immutable-polynomial.jl")
+include("standard-basis/sparse-polynomial.jl")
+
+include("promotions.jl")
+
 include("polynomials/ngcd.jl")
 include("polynomials/multroot.jl")
-include("polynomials/ChebyshevT.jl")
+
+include("polynomials/chebyshev.jl")
+
 
 # Rational functions
 include("rational-functions/common.jl")
@@ -31,7 +48,6 @@ include("rational-functions/fit.jl")
 #include("rational-functions/rational-transfer-function.jl")
 include("rational-functions/plot-recipes.jl")
 
-
 # compat; opt-in with `using Polynomials.PolyCompat`
 include("polynomials/Poly.jl")