From 65f8e077cdb16c4f152a3391da0cac76ac139ba7 Mon Sep 17 00:00:00 2001 From: jverzani Date: Fri, 28 Jul 2023 14:22:34 -0400 Subject: [PATCH] WIP: doc adjustments --- docs/src/extending.md | 187 ++++++++++++++++++++++- src/promotions.jl | 8 + src/show.jl | 5 +- src/standard-basis/standard-basis.jl | 6 +- src/standard-basis/standard-dense.jl | 95 ++++++++++++ src/standard-basis/standard-immutable.jl | 55 +++++++ src/standard-basis/standard-sparse.jl | 43 ++++++ 7 files changed, 394 insertions(+), 5 deletions(-) diff --git a/docs/src/extending.md b/docs/src/extending.md index 9a815044..e7e44264 100644 --- a/docs/src/extending.md +++ b/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,188 @@ 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 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{α}}}) where {α} = + "L^$(α)" +``` + +The basis symbol has no default. 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) +``` + +```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. diff --git a/src/promotions.jl b/src/promotions.jl index d43fbb3d..cb4b3c05 100644 --- a/src/promotions.jl +++ b/src/promotions.jl @@ -2,6 +2,14 @@ Base.promote_rule(::Type{P}, ::Type{Q}) where {B,T,S,X, P<:AbstractUnivariatePolynomial{B,T,X}, Q<:AbstractUnivariatePolynomial{B,S,X}} = MutableDensePolynomial{B,promote_type(T,S),X} +Base.promote_rule(::Type{P}, ::Type{Q}) where {B,T,S,X, + P<:AbstractPolynomial{T,X}, + Q<:AbstractUnivariatePolynomial{B,S,X}} = MutableDensePolynomial{B,promote_type(T,S),X} + +Base.promote_rule(::Type{P}, ::Type{Q}) where {B,T,S,X, + P<:AbstractUnivariatePolynomial{B,T,X}, + Q<:AbstractPolynomial{S,X}} = MutableDensePolynomial{B,promote_type(T,S),X} + # Methods to ensure that matrices of polynomials behave as desired Base.promote_rule(::Type{P},::Type{Q}) where {T,X, P<:AbstractPolynomial{T,X}, S, Q<:AbstractPolynomial{S,X}} = diff --git a/src/show.jl b/src/show.jl index 0f773ed4..61b7f338 100644 --- a/src/show.jl +++ b/src/show.jl @@ -59,10 +59,13 @@ showone(::Type{<:AbstractPolynomial{S}}) where {S} = false Common Printing =# +_typealias(::Type{P}) where {P<:AbstractPolynomial} = P.name.wrapper # allows for override + Base.show(io::IO, p::AbstractPolynomial) = show(io, MIME("text/plain"), p) function Base.show(io::IO, mimetype::MIME"text/plain", p::P) where {P<:AbstractPolynomial} - print(io,"$(P.name.wrapper)(") + print(io, _typealias(P)) + print(io, "(") printpoly(io, p, mimetype) print(io,")") end diff --git a/src/standard-basis/standard-basis.jl b/src/standard-basis/standard-basis.jl index 6f2db9f6..b23eb4ce 100644 --- a/src/standard-basis/standard-basis.jl +++ b/src/standard-basis/standard-basis.jl @@ -1,6 +1,6 @@ struct StandardBasis <: AbstractBasis end -basis_symbol(::Type{<:AbstractUnivariatePolynomial{StandardBasis}}) = "x" +basis_symbol(::Type{P}) where {P<:AbstractUnivariatePolynomial{StandardBasis}} = string(indeterminate(P)) function printbasis(io::IO, ::Type{P}, j::Int, m::MIME) where {B<:StandardBasis, P <: AbstractUnivariatePolynomial{B}} iszero(j) && return # no x^0 print(io, basis_symbol(P)) @@ -28,7 +28,7 @@ end function Base.convert(P::Type{PP}, q::Q) where {PP <: StandardBasisPolynomial, B<:StandardBasis,T,X, Q<:AbstractUnivariatePolynomial{B,T,X}} minimumexponent(P) > firstindex(chop(q)) && - throw(ArgumentError("Degree of polynomial less than minimum degree of polynomial type $(⟒(P))")) + throw(ArgumentError("Lowest degree term of polynomial less than the minimum degree of the polynomial type $(⟒(P))")) isa(q, PP) && return p T′ = _eltype(P,q) @@ -86,7 +86,7 @@ function integrate(p::AbstractUnivariatePolynomial{B,T,X}) where {B <: StandardB cs = _zeros(p, z, N+1) os = offset(p) @inbounds for (i, cᵢ) ∈ pairs(p) - i == -1 && (iszero(cᵢ) ? continue : throw(ArgumentError("Laurent polynomial with 1/x term"))) + i == -1 && (iszero(cᵢ) ? continue : throw(ArgumentError("Can't integrate Laurent polynomial with `x⁻¹` term"))) #cs[i + os] = cᵢ / (i+1) cs = _set(cs, i + 1 + os, cᵢ / (i+1)) end diff --git a/src/standard-basis/standard-dense.jl b/src/standard-basis/standard-dense.jl index b463cdc6..8ed0f28f 100644 --- a/src/standard-basis/standard-dense.jl +++ b/src/standard-basis/standard-dense.jl @@ -4,9 +4,104 @@ #const Polynomial = MutableDensePolynomial{StandardBasis} #export Polynomial +""" + LaurentPolynomial{T,X}(coeffs::AbstractVector, [m::Integer = 0], [var = :x]) + +A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of the form `a_{m}x^m + ... + a_{n}x^n` where `m,n` are integers (not necessarily positive) with ` m <= n`. + +The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. +The argument `m` represents the lowest exponent of the variable in the series, and is taken to be zero by default. + +Laurent polynomials and standard basis polynomials promote to Laurent polynomials. Laurent polynomials may be converted to a standard basis polynomial when `m >= 0` +. + +Integration will fail if there is a `x⁻¹` term in the polynomial. + +!!! note + `LaurentPolynomial` is an alias for `MutableDensePolynomial{StandardBasis}`. + +!!! note + `LaurentPolynomial` is axis-aware, unlike the other polynomial types in this package. + +# Examples: +```jldoctest laurent +julia> using Polynomials + +julia> P = LaurentPolynomial; + +julia> p = P([1,1,1], -1) +LaurentPolynomial(x⁻¹ + 1 + x) + +julia> q = P([1,1,1]) +LaurentPolynomial(1 + x + x²) + +julia> pp = Polynomial([1,1,1]) +Polynomial(1 + x + x^2) + +julia> p + q +LaurentPolynomial(x⁻¹ + 2 + 2*x + x²) + +julia> p * q +LaurentPolynomial(x⁻¹ + 2 + 3*x + 2*x² + x³) + +julia> p * pp +LaurentPolynomial(x⁻¹ + 2 + 3*x + 2*x² + x³) + +julia> pp - q +LaurentPolynomial(0) + +julia> derivative(p) +LaurentPolynomial(-x⁻² + 1) + +julia> integrate(q) +LaurentPolynomial(1.0*x + 0.5*x² + 0.3333333333333333*x³) + +julia> integrate(p) # x⁻¹ term is an issue +ERROR: ArgumentError: Can't integrate Laurent polynomial with `x⁻¹` term + +julia> integrate(P([1,1,1], -5)) +LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²) + +julia> x⁻¹ = inv(variable(LaurentPolynomial)) # `inv` defined on monomials +LaurentPolynomial(1.0*x⁻¹) + +julia> p = Polynomial([1,2,3]) +Polynomial(1 + 2*x + 3*x^2) + +julia> x = variable() +Polynomial(x) + +julia> x^degree(p) * p(x⁻¹) # reverses coefficients +LaurentPolynomial(3.0 + 2.0*x + 1.0*x²) +``` +""" const LaurentPolynomial = MutableDensePolynomial{StandardBasis} export LaurentPolynomial +_typealias(::Type{P}) where {P<:LaurentPolynomial} = "LaurentPolynomial" + +# how to show term. Only needed here to get unicode exponents to match old LaurentPolynomial.jl type +function showterm(io::IO, ::Type{P}, pj::T, var, j, first::Bool, mimetype) where {T, P<:MutableDensePolynomial{StandardBasis,T}} + if _iszero(pj) return false end + + pj = printsign(io, pj, first, mimetype) + if hasone(T) + if !(_isone(pj) && !(showone(T) || j == 0)) + printcoefficient(io, pj, j, mimetype) + end + else + printcoefficient(io, pj, j, mimetype) + end + + iszero(j) && return true + printproductsign(io, pj, j, mimetype) + print(io, indeterminate(P)) + j == 1 && return true + unicode_exponent(io, j) # print(io, exponent_text(j, mimetype)) + return true +end + + function evalpoly(c, p::MutableDensePolynomial{StandardBasis,T,X}) where {T,X} iszero(p) && return zero(T) * zero(c) EvalPoly.evalpoly(c, p.coeffs) * c^p.order[] diff --git a/src/standard-basis/standard-immutable.jl b/src/standard-basis/standard-immutable.jl index b61dd1d5..fdc5249c 100644 --- a/src/standard-basis/standard-immutable.jl +++ b/src/standard-basis/standard-immutable.jl @@ -1,7 +1,62 @@ ## Immutable dense / standard basis specific polynomial code +""" + ImmutablePolynomial{T, X, N}(coeffs) + +Construct an immutable (static) polynomial from its coefficients +`a₀, a₁, …, aₙ`, +lowest order first, optionally in terms of the given variable `x` +where `x` can be a character, symbol, or string. + +If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct +this through `ImmutablePolynomial((a_0, a_1, ..., a_n))` (assuming +`a_n ≠ 0`). As well, a vector or number can be used for construction. + + +The usual arithmetic operators are overloaded to work with polynomials +as well as with combinations of polynomials and scalars. However, +operations involving two non-constant polynomials of different variables causes an +error. Unlike other polynomials, `setindex!` is not defined for `ImmutablePolynomials`. + +As the degree of the polynomial (`+1`) is a compile-time constant, +several performance improvements are possible. For example, immutable +polynomials can take advantage of faster polynomial evaluation +provided by `evalpoly` from Julia 1.4; similar methods are also used +for addition and multiplication. + +However, as the degree is included in the type, promotion between +immutable polynomials can not promote to a common type. As such, they +are precluded from use in rational functions. + +!!! note + `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first + index always corresponding to the constant term. + +# Examples + +```jldoctest +julia> using Polynomials + +julia> ImmutablePolynomial((1, 0, 3, 4)) +ImmutablePolynomial(1 + 3*x^2 + 4*x^3) + +julia> ImmutablePolynomial((1, 2, 3), :s) +ImmutablePolynomial(1 + 2*s + 3*s^2) + +julia> one(ImmutablePolynomial) +ImmutablePolynomial(1.0) +``` + +!!! note + This was modeled after [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) by `@tkoolen`. + +!!! note + `ImmutablePolynomial` is an alias for `ImmutableDensePolynomial{StandardBasis}`. + +""" ImmutablePolynomial = ImmutableDensePolynomial{StandardBasis} export ImmutablePolynomial +_typealias(::Type{P}) where {P<:ImmutablePolynomial} = "ImmutablePolynomial" evalpoly(x, p::ImmutableDensePolynomial{B,T,X,0}) where {B<:StandardBasis,T,X} = zero(T)*zero(x) evalpoly(x, p::ImmutableDensePolynomial{B,T,X,N}) where {B<:StandardBasis,T,X,N} = EvalPoly.evalpoly(x, p.coeffs) diff --git a/src/standard-basis/standard-sparse.jl b/src/standard-basis/standard-sparse.jl index f9a02e7d..2a20ca13 100644 --- a/src/standard-basis/standard-sparse.jl +++ b/src/standard-basis/standard-sparse.jl @@ -1,7 +1,50 @@ ## Standard basis + sparse storage + +""" + SparsePolynomial{T, X}(coeffs::Dict{Int,T}) + +Polynomials in the standard basis backed by a dictionary holding the +non-zero coefficients. For polynomials of high degree, this might be +advantageous. + +# Examples: + +```jldoctest +julia> using Polynomials + +julia> P = SparsePolynomial; + +julia> p, q = P([1,2,3]), P([4,3,2,1]) +(SparsePolynomial(1 + 2*x + 3*x^2), SparsePolynomial(4 + 3*x + 2*x^2 + x^3)) + +julia> p + q +SparsePolynomial(5 + 5*x + 5*x^2 + x^3) + +julia> p * q +SparsePolynomial(4 + 11*x + 20*x^2 + 14*x^3 + 8*x^4 + 3*x^5) + +julia> p + 1 +SparsePolynomial(2 + 2*x + 3*x^2) + +julia> q * 2 +SparsePolynomial(8 + 6*x + 4*x^2 + 2*x^3) + +julia> p = Polynomials.basis(P, 10^9) - Polynomials.basis(P,0) # also P(Dict(0=>-1, 10^9=>1)) +SparsePolynomial(-1.0 + 1.0*x^1000000000) + +julia> p(1) +0.0 +``` + +!!! note + `SparsePolynomial` is an alias for `MutableSparsePolynomial{StandardBasis}`. + +""" const SparsePolynomial = MutableSparsePolynomial{StandardBasis} # const is important! export SparsePolynomial +_typealias(::Type{P}) where {P<:SparsePolynomial} = "SparsePolynomial" + function evalpoly(x, p::MutableSparsePolynomial) tot = zero(p[0]*x)