V2.0.0 (#310)

* move symbol into type

* change what iteration means for a polynomial

* adjust integrate to return polynomial of same type, even if NaN

* clean up offset vector tests

* drop OffsetArrays dependency; clean up scalar operations for ImmutablePolynomial

* change some errors into ArgumentError

* replace isnothing; standardize check for same variable

* registerN remove, previously deprecated

* update documentation, add example to `extending.md`

* clean up integration, base on `integrate(p)`, not `integrate(p,c)`

* work on promotion for + and * operations; move computations into `common` for different storage types

* laurent polynomial depwarns removed

* fix basis symbol handling

* adjust evaluation so that implemented evalpoly instead of call syntax is expected.

* add example to extending; specialize evalpoly rather than call syntax

* refactor truncate! chop! move computations based on storage type into `common.jl`

* move things into common (zero, one,variable, basis)

* fix LaurentPolynomial evaluation to use faster `evalpoly`

* regularize treatment of indeterminate in zero,one,variable,basis

Author: jverzani

Project.toml

@@ -2,25 +2,25 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "1.2.0"
+version = "2.0.0"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 
+
 [compat]
 Intervals = "0.5, 1.0, 1.3"
-OffsetArrays = "1"
 RecipesBase = "0.7, 0.8, 1"
 julia = "1"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "SparseArrays", "SpecialFunctions", "Test"]
+test = ["LinearAlgebra", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]

README.md

@@ -208,7 +208,9 @@ Polynomial objects also have other methods:
 
 * [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
 
-* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) and [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
+* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl), [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series, [Hecke.jl](https://github.com/thofma/Hecke.jl) for algebraic number theory.
+
+* [CommutativeAlgebra](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
 
 * [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 

docs/src/extending.md

@@ -21,14 +21,129 @@ As always, if the default implementation does not work or there are more efficie
 | Constructor | x | |
 | Type function (`(::P)(x)`) | x | |
 | `convert(::Polynomial, ...)` | | Not required, but the library is built off the [`Polynomial`](@ref) type, so all operations are guaranteed to work with it. Also consider writing the inverse conversion method. |
+| `Polynomials.evalpoly(x, p::P)` |  to evaluate the polynomial at `x` (`Base.evalpoly` okay post `v"1.4.0"`) |
 | `domain` | x | Should return an  [`AbstractInterval`](https://invenia.github.io/Intervals.jl/stable/#Intervals-1) |
 | `vander` | | Required for [`fit`](@ref) |
 | `companion` | | Required for [`roots`](@ref) |
-| `fromroots` | | By default, will form polynomials using `prod(variable(::P) - r)` for reach root `r`|
-| `+(::P, ::P)` | | Addition of polynomials |
-| `-(::P, ::P)` | | Subtraction of polynomials |
 | `*(::P, ::P)` | | Multiplication of polynomials |
 | `divrem` | | Required for [`gcd`](@ref)|
+| `one`| | Convenience to find constant in new basis |
 | `variable`| | Convenience to find monomial `x` in new  basis|
 
-Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended.
+Check out both the [`Polynomial`](@ref) and [`ChebyshevT`](@ref) for examples of this interface being extended. 
+
+## Example
+
+The following shows a minimal example where the polynomial aliases the vector defining the coefficients. 
+The constructor ensures that there are no trailing zeros. The `@register` call ensures a common interface. This example subtypes `StandardBasisPolynomial`, not `AbstractPolynomial`, and consequently inherits the methods above that otherwise would have been required. For other bases,  more methods may be necessary to define  (again, refer to [`ChebyshevT`](@ref) for an example).
+
+```jldoctest AliasPolynomial
+julia> using Polynomials
+
+julia> struct AliasPolynomial{T <: Number, X} <: Polynomials.StandardBasisPolynomial{T, X}
+                  coeffs::Vector{T}
+                  function AliasPolynomial{T, X}(coeffs::Vector{S}) where {T, X, S}
+                      p = new{T,X}(coeffs)
+                      chop!(p)
+                  end
+              end
+
+julia> Polynomials.@register AliasPolynomial
+```
+
+To see this new polynomial type in action, we have:
+
+```jldoctest AliasPolynomial
+julia> xs = [1,2,3,4];
+
+julia> p = AliasPolynomial(xs)
+AliasPolynomial(1 + 2*x + 3*x^2 + 4*x^3)
+
+julia> q = AliasPolynomial(1.0, :y)
+AliasPolynomial(1.0)
+
+julia> p + q
+AliasPolynomial(2.0 + 2.0*x + 3.0*x^2 + 4.0*x^3)
+
+julia> p * p
+AliasPolynomial(1 + 4*x + 10*x^2 + 20*x^3 + 25*x^4 + 24*x^5 + 16*x^6)
+
+julia> (derivative ∘ integrate)(p) == p
+true
+
+julia> p(3)
+142
+```
+
+For the `Polynomial` type, the default on operations is to copy the array. For this type, it might seem reasonable -- to avoid allocations -- to update the coefficients in place for scalar addition and scalar multiplication. 
+
+Scalar addition, `p+c`, defaults to `p + c*one(p)`, or polynomial addition, which is not inplace without addition work. As such, we create a new method and an infix operator
+
+```jldoctest AliasPolynomial
+julia> function scalar_add!(p::AliasPolynomial{T}, c::T) where {T}
+           p.coeffs[1] += c
+           p
+       end;
+
+julia> p::AliasPolynomial ⊕ c::Number = scalar_add!(p,c);
+
+```
+
+Then we have:
+
+```jldoctest AliasPolynomial
+julia> p
+AliasPolynomial(1 + 2*x + 3*x^2 + 4*x^3)
+
+julia> p ⊕ 2
+AliasPolynomial(3 + 2*x + 3*x^2 + 4*x^3)
+
+julia> p
+AliasPolynomial(3 + 2*x + 3*x^2 + 4*x^3)
+```
+
+The viewpoint that a polynomial represents a vector of coefficients  leads to an expectation that vector operations should match when possible. Scalar multiplication is a vector operation, so it seems reasonable to override the broadcast machinery to implement an in place operation (e.g. `p .*= 2`). By default, the polynomial types are not broadcastable over their coefficients. We would need to make a change there and modify the `copyto!` function:
+
+
+```jldoctest AliasPolynomial
+julia> Base.broadcastable(p::AliasPolynomial) = p.coeffs;
+
+
+julia> Base.ndims(::Type{<:AliasPolynomial}) = 1
+
+
+julia> Base.copyto!(p::AliasPolynomial, x) = (copyto!(p.coeffs, x); chop!(p));
+
+```
+
+The last `chop!` call would ensure that there are no trailing zeros in the coefficient vector after multiplication, as multiplication by `0` is possible.
+
+Then we might have:
+
+```jldoctest AliasPolynomial
+julia> p
+AliasPolynomial(3 + 2*x + 3*x^2 + 4*x^3)
+
+julia> p .*= 2
+AliasPolynomial(6 + 4*x + 6*x^2 + 8*x^3)
+
+julia> p
+AliasPolynomial(6 + 4*x + 6*x^2 + 8*x^3)
+
+julia> p ./= 2
+AliasPolynomial(3 + 2*x + 3*x^2 + 4*x^3)
+```
+
+Trying to divide again would throw an error, as the result would not fit with the integer type of `p`. 
+
+Now `p` is treated as the vector `p.coeffs`, as regards broadcasting, so some things may be surprising, for example this expression returns a vector, not a polynomial:
+
+```jldoctest AliasPolynomial
+julia> p .+ 2
+4-element Array{Int64,1}:
+ 5
+ 4
+ 5
+ 6
+```
+

docs/src/index.md

@@ -97,7 +97,7 @@ julia> q = Polynomial([1, 2, 3], :s)
 Polynomial(1 + 2*s + 3*s^2)
 
 julia> p + q
-ERROR: Polynomials must have same variable
+ERROR: ArgumentError: Polynomials have different indeterminates
 [...]
 ```
 
@@ -199,7 +199,7 @@ julia> derivative(p1)
 ChebyshevT(2.0⋅T_0(x) + 12.0⋅T_1(x))
 
 julia> integrate(p2)
-ChebyshevT(0.25⋅T_0(x) - 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))
+ChebyshevT(- 1.0⋅T_1(x) + 0.25⋅T_2(x) + 0.3333333333333333⋅T_3(x))
 
 julia> convert(Polynomial, p1)
 Polynomial(-2.0 + 2.0*x + 6.0*x^2)
@@ -213,15 +213,68 @@ ChebyshevT(2.5⋅T_0(x) + 2.0⋅T_1(x) + 1.5⋅T_2(x))
 
 ### Iteration
 
-If its basis is implicit, then a polynomial may be  seen as just a vector of  coefficients. Vectors or 1-based, but, for convenience, polynomial types are 0-based, for purposes of indexing (e.g. `getindex`, `setindex!`, `eachindex`). Iteration over a polynomial steps through the basis vectors, e.g. `a_0`, `a_1*x`, ...
+If its basis is implicit, then a polynomial may be  seen as just a vector of  coefficients. Vectors or 1-based, but, for convenience, polynomial types are 0-based, for purposes of indexing (e.g. `getindex`, `setindex!`, `eachindex`). Iteration over a polynomial steps through the underlying coefficients.
 
 ```jldoctest
 julia> as = [1,2,3,4,5]; p = Polynomial(as);
 
 julia> as[3], p[2], collect(p)[3]
-(3, 3, Polynomial(3*x^2))
+(3, 3, 3)
 ```
 
+
+The `pairs` iterator, iterates over the indices and coefficients, attempting to match how `pairs` applies to the underlying storage model:
+
+```jldoctest
+julia> v = [1,2,0,4]
+4-element Array{Int64,1}:
+ 1
+ 2
+ 0
+ 4
+
+julia> p,ip,sp,lp = Polynomial(v), ImmutablePolynomial(v), SparsePolynomial(v), LaurentPolynomial(v, -1);
+
+julia> collect(pairs(p))
+4-element Array{Pair{Int64,Int64},1}:
+ 0 => 1
+ 1 => 2
+ 2 => 0
+ 3 => 4
+
+julia> collect(pairs(ip)) == collect(pairs(p))
+true
+
+julia> collect(pairs(sp)) # unordered dictionary with only non-zero terms
+3-element Array{Pair{Int64,Int64},1}:
+ 0 => 1
+ 3 => 4
+ 1 => 2
+
+julia> collect(pairs(lp))
+4-element Array{Pair{Int64,Int64},1}:
+ -1 => 1
+  0 => 2
+  1 => 0
+  2 => 4
+```
+
+
+The unexported `monomials` iterator iterates over the terms (`p[i]*Polynomials.basis(p,i)`) of the polynomial:
+
+```jldoctest
+julia> p = Polynomial([1,2,0,4], :u)
+Polynomial(1 + 2*u + 4*u^3)
+
+julia> collect(Polynomials.monomials(p))
+4-element Array{Any,1}:
+ Polynomial(1)
+ Polynomial(2*u)
+ Polynomial(0)
+ Polynomial(4*u^3)
+```
+
+
 ## Related Packages
 
 * [StaticUnivariatePolynomials.jl](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Fixed-size univariate polynomials backed by a Tuple
@@ -234,12 +287,15 @@ julia> as[3], p[2], collect(p)[3]
 
 * [PolynomialRings](https://github.com/tkluck/PolynomialRings.jl) A library for arithmetic and algebra with multi-variable polynomials.
 
-* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl) and [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series
+* [AbstractAlgebra.jl](https://github.com/wbhart/AbstractAlgebra.jl), [Nemo.jl](https://github.com/wbhart/Nemo.jl) for generic polynomial rings, matrix spaces, fraction fields, residue rings, power series, [Hecke.jl](https://github.com/thofma/Hecke.jl) for algebraic number theory.
+
+* [CommutativeAlgebra](https://github.com/KlausC/CommutativeRings.jl) the start of a computer algebra system specialized to discrete calculations with support for polynomials.
 
 * [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) for a fast complex polynomial root finder. For larger degree problems, also [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) and [AMRVW](https://github.com/jverzani/AMRVW.jl).
 
 
 
+
 ## Contributing
 
 If you are interested in this project, feel free to open an issue or pull request! In general, any changes must be thoroughly tested, allow deprecation, and not deviate too far from the common interface. All PR's must have an updated project version, as well, to keep the continuous delivery cycle up-to-date.

docs/src/polynomials/chebyshev.md

@@ -31,7 +31,7 @@ ChebyshevT(1⋅T_0(x) + 3⋅T_2(x) + 4⋅T_3(x))
 
 
 julia> p = convert(Polynomial, c)
-Polynomial(-2 - 12*x + 6*x^2 + 16*x^3)
+Polynomial(-2.0 - 12.0*x + 6.0*x^2 + 16.0*x^3)
 
 julia> convert(ChebyshevT, p)
 ChebyshevT(1.0⋅T_0(x) + 3.0⋅T_2(x) + 4.0⋅T_3(x))

src/Polynomials.jl

@@ -3,7 +3,10 @@ module Polynomials
 #  using GenericLinearAlgebra ## remove for now. cf: https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/71#issuecomment-743928205
 using LinearAlgebra
 using Intervals
-using OffsetArrays
+
+if VERSION >= v"1.4.0"
+    import Base: evalpoly
+end
 
 include("abstract.jl")
 include("show.jl")