Several updates

README.md

@@ -5,47 +5,35 @@
 | [![][docs-stable-img]][docs-stable-url] | [![][pkg-0.5-img]][pkg-0.5-url] | [![Build Status][build-img]][build-url] [![Build Status][winbuild-img]][winbuild-url] |
 | [![][docs-latest-img]][docs-latest-url] | [![][pkg-0.6-img]][pkg-0.6-url] | [![Coveralls branch][coveralls-img]][coveralls-url] [![Codecov branch][codecov-img]][codecov-url] |
 
-Basic arithmetic, integration, differentiation and evaluation over sparse multivariate polynomials and sparse multivariate moments.
-Both commutative and non-commutative variables are supported.
-The following types are defined:
+This package provides an interface for manipulating multivariate polynomials.
+Implementing algorithms on polynomials using this interface will allow the algorithm to work for all polynomials implementing the interface.
 
-* `PolyVar{C}`: A variable which is commutative with `*` when `C` is `true`. Commutative variables are created using the `@polyvar` macro, e.g. `@polyvar x y`, `@polyvar x[1:8]` and non-commutative variables are created likewise using the `@ncpolyvar` macro.
-* `Monomial{C}`: A product of variables: e.g. `x*y^2`.
-* `Term{C, T}`: A product between an element of type `T` and a `Monomial{C}`, e.g `2x`, `3.0x*y^2`.
-* `Polynomial{C, T}`: A sum of `Term{C, T}`, e.g. `2x + 3.0x*y^2 + y`.
-* `Moment{C, T}`: The multivariate moment of type `T` of a measure, e.g. `E_μ[x*y^2]` is the moment of `μ` corresponding to the monomial `x*y^2`.
-* `Measure{C, T}`: A combination of `Moment{C, T}` of a measure, e.g. the moments of `x`, `x*y^2` and `y`.
+The interface contains functions for accessing the coefficients, monomials, terms of the polynomial, defines arithmetic operations on them, rational functions, division with remainder, calculus/differentiation and evaluation/substitution.
 
-All common algebraic operations between those types are designed to be as efficient as possible without doing any assumption on `T`.
-Typically, one imagine `T` to be a subtype of `Number` but it can be anything.
-This is useful for example in the package [PolyJuMP](https://github.com/JuliaOpt/PolyJuMP.jl) where `T` is often an affine expression of [JuMP](https://github.com/JuliaOpt/JuMP.jl) decision variables.
-The commutativity of `T` with `*` is not assumed, even if it is the coefficient of a monomial of commutative variables.
-However, commutativity of `T` and of the variables `+` is always assumed.
-This allows to keep the terms and moments always sorted (Graded Lexicographic order is used) in polynomial and measure which enables more efficient operations.
+The following packages provides multivariate polynomials that implement the interface.
+
+* [TypedPolynomials](https://github.com/rdeits/TypedPolynomials.jl) : Commutative polynomials of arbitrary coefficient types
+* [DynamicPolynomials](https://github.com/blegat/DynamicPolynomials.jl) : Commutative and non-commutative polynomials of arbitrary coefficient types
 
 Below is a simple usage example
 ```julia
-@polyvar x y
+using TypedPolynomials
+@polyvar x y # assigns x (resp. y) to a variable of name x (resp. y)
 p = 2x + 3.0x*y^2 + y
-differentiate(p, x) # compute the derivative of p with respect to x
-differentiate(p, [x, y]) # compute the gradient of p
-p([y, x], [x, y]) # replace any x by y and y by x
-subs(p, [x^2], [y]) # replace any occurence of y by x^2
-p([1, 2], [x, y]) # evaluate p at [1, 2]
+@test differentiate(p, x) # compute the derivative of p with respect to x
+@test differentiate.(p, (x, y)) # compute the gradient of p
+@test p((x, y)=>(y, x)) # replace any x by y and y by x
+@test subs(p, y=>x^2) # replace any occurence of y by x^2
+@test p(x=>1, y=>2) # evaluate p at [1, 2]
 ```
-Below is an example with `@polyvar x[1:n]`
+Below is an example with `@polyvar x[1:3]`
 ```julia
-n = 3
 A = rand(3, 3)
-@polyvar x[1:n]
-p = dot(x, x) # x_1^2 + x_2^2 + x_3^2
-p(A*x, x) # corresponds to dot(A*x, A*x)
-subs(p, [2, 3], [x[1], x[3]]) # x_2^2 + 13
+@polyvar x[1:3] # assign x to a tuple of variables x1, x2, x3
+p = sum(x .* x) # x_1^2 + x_2^2 + x_3^2
+subs(p, x[1]=>2, x[3]=>3) # x_2^2 + 13
+p(x=>A*vec(x)) # corresponds to dot(A*x, A*x), need vec to convert the tuple to a vector
 ```
-Note that, when doing substitution, it is required to give the `PolyVar` ordering that is meant.
-Indeed, the ordering between the `PolyVar` is not alphabetical but rather by order of creation
-which can be undeterministic with parallel computing.
-Therefore, this order cannot be used for substitution, even as a default (see [here](https://github.com/blegat/MultivariatePolynomials.jl/issues/3) for a discussion about this).
 
 ## See also
 

src/div.jl

@@ -12,7 +12,11 @@ function divides(t1::AbstractTermLike, t2::AbstractTermLike)
 end
 divides(t1::AbstractVariable, t2::AbstractVariable) = t1 == t2
 
+Base.gcd(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(min, m1, m2)
+Base.lcm(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(max, m1, m2)
+
 # _div(a, b) assumes that b divides a
+_div(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(-, m1, m2)
 function _div(t::AbstractTerm, m::AbstractMonomial)
     coefficient(t) * _div(monomial(t), m)
 end

src/mono.jl

@@ -1,4 +1,6 @@
-export name, constantmonomial, emptymonovec, monovec, monovectype, sortmonovec, mergemonovec
+export name, constantmonomial, emptymonovec, monovec, monovectype, sortmonovec, mergemonovec, mapexponents
+
+mapexponents(f, m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(f, monomial(m1), monomial(m2))
 
 Base.copy(x::AbstractVariable) = x
 

src/operators.jl

@@ -47,7 +47,7 @@ multconstant(p::AbstractPolynomialLike, α) = multconstant(polynomial(p), α)
 multconstant(α, t::AbstractTermLike)    = (α*coefficient(t)) * monomial(t)
 multconstant(t::AbstractTermLike, α)    = (coefficient(t)*α) * monomial(t)
 
-(*)(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(promote(m1, m2)...)
+(*)(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(+, m1, m2)
 #(*)(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(monomial(m1), monomial(m2))
 
 (*)(m::AbstractMonomialLike, t::AbstractTermLike) = coefficient(t) * (m * monomial(t))

src/substitution.jl

@@ -51,7 +51,7 @@ _polynomial(α) = α
 _polynomial(p::APL) = polynomial(p)
 function substitute(st::AST, p::AbstractPolynomial, s::Substitutions)
     if iszero(p)
-        _polynomial(substitute(st, zero(termtype(p)), s))
+        _polynomial(substitute(st, zeroterm(p), s))
     else
         ts = terms(p)
         r1 = substitute(st, ts[1], s)

src/term.jl

@@ -10,9 +10,17 @@ function Base.hash(t::AbstractTerm, u::UInt)
     end
 end
 
-Base.zero(::Type{TT}) where {T, TT<:AbstractTermLike{T}} = zero(T) * constantmonomial(TT)
+# zero should return a polynomial since it is often used to keep the result of a summation of terms.
+# For example, Base.vecdot(x::Vector{<:AbstractTerm}, y:Vector{Int}) starts with `s = zero(dot(first(x), first(y)))` and then adds terms.
+# We want `s` to start as a polynomial for this operation to be efficient.
+#Base.zero(::Type{TT}) where {T, TT<:AbstractTermLike{T}} = zero(T) * constantmonomial(TT)
+#Base.zero(t::AbstractTermLike{T}) where {T} = zero(T) * constantmonomial(t)
+zeroterm(::Type{TT}) where {T, TT<:APL{T}} = zero(T) * constantmonomial(TT)
+zeroterm(t::APL{T}) where {T} = zero(T) * constantmonomial(t)
+
+Base.zero(::Type{TT}) where {T, TT<:AbstractTermLike{T}} = zero(polynomialtype(TT))
+Base.zero(t::AbstractTermLike{T}) where {T} = zero(polynomialtype(t))
 Base.one(::Type{TT}) where {T, TT<:AbstractTermLike{T}} = one(T) * constantmonomial(TT)
-Base.zero(t::AbstractTermLike{T}) where {T} = zero(T) * constantmonomial(t)
 Base.one(t::AbstractTermLike{T}) where {T} = one(T) * constantmonomial(t)
 
 monomial(m::AbstractMonomial) = m
@@ -36,9 +44,9 @@ When applied on a polynomial, it throws an error if it has more than one term.
 When applied to a term, it is the identity and does not copy it.
 When applied to a monomial, it create a term of type `AbstractTerm{Int}`.
 """
-function term(p::APL)
+function term(p::APL{T}) where T
     if nterms(p) == 0
-        zero(termtype(p))
+        zeroterm(p)
     elseif nterms(p) > 1
         error("A polynomial is more than one term cannot be converted to a term.")
     else

test/div.jl

@@ -1,4 +1,13 @@
 @testset "Division" begin
+    @testset "GCD and LCM" begin
+        Mod.@polyvar x y z
+        @test gcd(x*y, y) == y
+        @test lcm(x*y, y) == x*y
+        @test gcd(x*y^2*z, y*z^3) == y*z
+        @test lcm(x*y^2*z, y*z^3) == x*y^2*z^3
+        @test gcd(x^2*y^7*z^3, x^4*y^5*z^2) == x^2*y^5*z^2
+        @test lcm(x^2*y^7*z^3, x^4*y^5*z^2) == x^4*y^7*z^3
+    end
     # Taken from
     # Ideals, Varieties, and Algorithms
     # Cox, Little and O'Shea, Fourth edition

test/mono.jl

@@ -17,7 +17,7 @@ const MP = MultivariatePolynomials
         @test !iszero(x)
         @test zero(x) == 0
         @test iszero(zero(x))
-        @test zero(x) isa AbstractTerm{Int}
+        @test zero(x) isa AbstractPolynomial{Int}
         @inferred zero(x)
         @test one(x) == 1
         @test one(x) isa AbstractMonomial
@@ -36,7 +36,7 @@ const MP = MultivariatePolynomials
     @testset "Monomial" begin
         Mod.@polyvar x
         @test zero(x^2) == 0
-        @test zero(x^2) isa AbstractTerm{Int}
+        @test zero(x^2) isa AbstractPolynomial{Int}
         @inferred zero(x^2)
         @test one(x^2) == 1
         @test one(x^2) isa AbstractMonomial