WIP: doc adjustments

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.

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}} =

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

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

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[]