Factored polynomial (#340)

* add factored polynomial type

* update README, docs, fix bug in divrem

* work around v1.0-1.2

Author: jverzani

README.md

@@ -19,7 +19,9 @@ Basic arithmetic, integration, differentiation, evaluation, and root finding ove
 * `ImmutablePolynomial` –⁠ standard basis polynomials backed by a [Tuple type](https://docs.julialang.org/en/v1/manual/functions/#Tuples-1) for faster evaluation of values
 * `SparsePolynomial` –⁠ standard basis polynomial backed by a [dictionary](https://docs.julialang.org/en/v1/base/collections/#Dictionaries-1) to hold  sparse high-degree  polynomials
 * `LaurentPolynomial` –⁠ [Laurent polynomials](https://docs.julialang.org/en/v1/base/collections/#Dictionaries-1), `a(x) = aₘ xᵐ + … + aₙ xⁿ` `m ≤ n`, `m,n ∈ ℤ` backed by an [offset array](https://github.com/JuliaArrays/OffsetArrays.jl); for example, if `m<0` and `n>0`, `a(x) = aₘ xᵐ + … + a₋₁ x⁻¹ + a₀ + a₁ x + … +  aₙ xⁿ`
+* `FactoredPolynomial` –⁠ standard basis polynomials, storing the roots, with multiplicity, and leading coefficient of a polynomial
 * `ChebyshevT` –⁠ [Chebyshev polynomials](https://en.wikipedia.org/wiki/Chebyshev_polynomials) of the first kind
+* `RationalFunction` - a type for ratios of polynomials.
 
 ## Usage
 

docs/src/index.md

@@ -159,6 +159,19 @@ julia> roots(Polynomial([0, 0, 1]))
  0.0
 ```
 
+For polynomials with multplicities, the non-exported `Polynomials.Multroot.multroot` method can avoid some numerical issues that `roots` will have. 
+
+The `FactoredPolynomial` type stores the roots (with multiplicities) and the leading coefficent of a polynomial. In this example, the `multroot` method is used internally to identify the roots of `p` below, in the conversion from the `Polynomial` type to the `FactoredPolynomial` type:
+
+```jldoctest
+julia> p = Polynomial([24, -50, 35, -10, 1])
+Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)
+
+julia> q = convert(FactoredPolynomial, p) # noisy form of `factor`:
+FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))
+```
+
+
 ### Fitting arbitrary data
 
 Fit a polynomial (of degree `deg`) to `x` and `y` using polynomial interpolation or a (weighted) least-squares approximation.
@@ -274,6 +287,25 @@ julia> collect(Polynomials.monomials(p))
  Polynomial(4*u^3)
 ```
 
+The `map` function for polynomials is idiosyncratic, as iteration over
+polynomials is over the vector of coefficients, but `map` will also
+maintain the type of the polynomial. Here we use `map` to smooth out
+the round-off error coming from the root-finding algorithm used
+internally when converting to the `FactoredPolynomial` type:
+
+
+```jldoctest
+julia> p = Polynomial([24, -50, 35, -10, 1])
+Polynomial(24 - 50*x + 35*x^2 - 10*x^3 + x^4)
+
+julia> q = convert(FactoredPolynomial, p) # noisy form of `factor`:
+FactoredPolynomial((x - 4.0000000000000036) * (x - 2.9999999999999942) * (x - 1.0000000000000002) * (x - 2.0000000000000018))
+
+julia> map(round, q, digits=10)
+FactoredPolynomial((x - 4.0) * (x - 2.0) * (x - 3.0) * (x - 1.0))
+```
+
+
 ## Relationship between the `T` and `P{T,X}`
 
 The addition of a polynomial and a scalar, such as
@@ -344,7 +376,7 @@ ERROR: ArgumentError: Polynomials have different indeterminates
 [...]
 ```
 
-an error thrown.
+an error is thrown.
 
 In general, arrays with mixtures of non-constant polynomials with *different* indeterminates will error. By default, an error will occur when constant polynomials with different indeterminates are used as components. However, for *typed* arrays, conversion will allow such constructs to be used.
 
@@ -388,6 +420,53 @@ julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}}
 Though were a non-constant polynomial with indeterminate `y` replacing
 `2one(q)` above, that addition would throw an error.
 
+## Rational functions
+
+The package provides support for rational functions -- fractions of polynomials (for most types). The construction of the basic type mirrors the construction of rational numbers.
+
+```jldoctest
+julia> P = FactoredPolynomial
+FactoredPolynomial
+
+julia> p,q = fromroots(P, [1,2,3,4]), fromroots(P, [2,2,3,5])
+(FactoredPolynomial((x - 4) * (x - 2) * (x - 3) * (x - 1)), FactoredPolynomial((x - 5) * (x - 2)² * (x - 3)))
+
+julia> pq = p // q
+((x - 4) * (x - 2) * (x - 3) * (x - 1)) // ((x - 5) * (x - 2)² * (x - 3))
+
+julia> lowest_terms(pq)
+((x - 4.0) * (x - 1.0)) // ((x - 5.0) * (x - 2.0))
+
+julia> d,r = residues(pq); r
+Dict{Float64, Vector{Float64}} with 2 entries:
+  5.0 => [1.33333]
+  2.0 => [0.666667]
+
+julia> x = variable(p);
+
+julia> for (λ, rs) ∈ r # reconstruct p/q from output of `residues`
+           for (i,rᵢ) ∈ enumerate(rs)
+               d += rᵢ//(x-λ)^i
+           end
+       end
+
+julia> d
+((x - 4.0) * (x - 1.0000000000000002)) // ((x - 5.0) * (x - 2.0))
+```
+
+A basic plot recipe is provided.
+
+```@example
+using Plots, Polynomials
+P = FactoredPolynomial
+p,q = fromroots(P, [1,2,3]), fromroots(P, [2,3,3,0])
+plot(p//q)
+savefig("rational_function.svg"); nothing # hide
+```
+
+![](rational_function.svg)
+
+
 ## Related Packages
 
 * [StaticUnivariatePolynomials.jl](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) Fixed-size univariate polynomials backed by a Tuple

src/Polynomials.jl

@@ -25,6 +25,7 @@ include("polynomials/ImmutablePolynomial.jl")
 include("polynomials/SparsePolynomial.jl")
 include("polynomials/LaurentPolynomial.jl")
 include("polynomials/pi_n_polynomial.jl")
+include("polynomials/factored_polynomial.jl")
 include("polynomials/ngcd.jl")
 include("polynomials/multroot.jl")
 

src/common.jl

@@ -996,6 +996,14 @@ function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{T};
     return r₀
 end
 
+"""
+    uvw(p,q; kwargs...)
+
+return `u` the gcd of `p` and `q`, and `v` and `w`, where `u*v = p` and `u*w = q`.
+"""
+uvw(p::AbstractPolynomial, q::AbstractPolynomial; kwargs...) = uvw(promote(p,q)...; kwargs...)
+uvw(p1::P, p2::P; kwargs...) where {P <:AbstractPolynomial} = throw(ArgumentError("uvw not defind"))
+
 """
     div(::AbstractPolynomial, ::AbstractPolynomial)
 """
@@ -1015,8 +1023,8 @@ Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 
 function Base.isapprox(p1::AbstractPolynomial{T,X},
-    p2::AbstractPolynomial{S,Y};
-    rtol::Real = (Base.rtoldefault(T, S, 0)),
+                       p2::AbstractPolynomial{S,Y};
+                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                        atol::Real = 0,) where {T,X,S,Y}
     assert_same_variable(p1, p2)
     (hasnan(p1) || hasnan(p2)) && return false  # NaN poisons comparisons