Fix #51

JuliaAlgebra · Aug 30, 2017 · 265b40e · 265b40e
1 parent 436714d
commit 265b40e
Show file tree

Hide file tree

Showing 22 changed files with 230 additions and 58 deletions.
diff --git a/README.md b/README.md
@@ -18,6 +18,7 @@ The following packages provides multivariate polynomials that implement the inte
 The following packages extend/use the interface:
 
 * [SemialgebraicSets](https://github.com/blegat/SemialgebraicSets.jl) : Sets defined by inequalities and equalities between polynomials and algorithms for solving polynomial systems of equations.
+* [HomotopyContinuation](https://github.com/saschatimme/HomotopyContinuation.jl) : Solving systems of polynomials via homotopy continuation.
 * [MultivariateMoments](https://github.com/blegat/MultivariateMoments.jl) : Moments of multivariate measures and their scalar product with polynomials.
 * [PolyJuMP](https://github.com/JuliaOpt/PolyJuMP.jl) : A [JuMP](https://github.com/JuliaOpt/JuMP.jl) extension for Polynomial Optimization.
 * [SumOfSquares](https://github.com/JuliaOpt/SumOfSquares.jl) : Certifying the nonnegativity of polynomials, minimizing/maximizing polynomials and optimization over sum of squares polynomials using Sum of Squares Programming.

diff --git a/docs/make.jl b/docs/make.jl
@@ -5,7 +5,10 @@ makedocs(
     sitename = "MultivariatePolynomials",
     pages = [
         "Introduction" => "index.md",
-        "Reference" => "apireference.md",
+        "Types" => "types.md",
+        "Substitution" => "substitution.md",
+        "Differentiation" => "differentiation.md",
+        "Division" => "division.md",
     ]
 )
 

diff --git a/docs/src/differentiation.md b/docs/src/differentiation.md
@@ -0,0 +1,8 @@
+# Differentiation
+
+Given a polynomial, say ``p(x, y) = 3x^2y + x + 2y + 1``, we can differentiate it by a variable, say ``x`` and get ``\partial p(x, y) / \partial x = 6xy + 1``.
+We can also differentiate it by both of its variable and get the vector ``[6xy+1, 3x^2+1]``.
+
+```@docs
+differentiate
+```
diff --git a/docs/src/division.md b/docs/src/division.md
@@ -0,0 +1,15 @@
+# Division
+
+The `gcd` and `lcm` functions of `Base` have been implemented for monomials, you have for example `gcd(x^2*y^7*z^3, x^4*y^5*z^2)` returning `x^2*y^5*z^2` and `lcm(x^2*y^7*z^3, x^4*y^5*z^2)` returning `x^4*y^7*z^3`.
+
+Given two polynomials, ``p`` and ``d``, there are unique ``r`` and ``q`` such that ``p = q d + r`` and the leading term of ``d`` does not divide the leading term of ``r``.
+You can obtain ``q`` using the `div` function and ``r`` using the `rem` function.
+The `divrem` function returns ``(q, r)``.
+
+Given a polynomial ``p`` and divisors ``d_1, \ldots, d_n``, one can find ``r`` and ``q_1, \ldots, q_n`` such that ``p = q_1 d_1 + \cdots + q_n d_n + r`` and none of the leading terms of ``q_1, \ldots, q_n`` divide the leading term of ``r``.
+You can obtain the vector ``[q_1, \ldots, q_n]`` using `div(p, d)` where ``d = [d_1, \ldots, d_n]`` and ``r`` using the `rem` function with the same arguments.
+The `divrem` function returns ``(q, r)``.
+
+```@docs
+divides
+```
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -1,18 +1,28 @@
 # MultivariatePolynomials
 
 [MultivariatePolynomials.jl](https://github.com/blegat/MultivariatePolynomials.jl) is an implementation independent library for manipulating multivariate polynomials.
-It defines abstract types and an API for multivariate monomials, terms, polynomials, moments and measures and gives default implementation for common operations on them using the API.
-If you want to manipulate multivariate polynomials easily and efficiently while being able to easily switch between different implementations, this library is exactly what you are looking for.
+It defines abstract types and an API for multivariate monomials, terms, polynomials and gives default implementation for common operations on them using the API.
 
-Supported operations are : basic arithmetic, rational polynomials, differentiation and evaluation/substitution, division and duality operations between polynomials and moments.
-There is also support for solving systems of equations (soon!) and building (semi)algebraic sets.
+On the one hand, This packages allows you to implement algorithms on multivariate polynomials that will be independant on the representation of the polynomial that will be chosen by the user.
+On the other hand, it allows the user to easily switch between different representations of polynomials to see which one is faster for the algorithm that he is using.
 
-Currently, the following implementations are available:
+Supported operations are : basic arithmetic, rational polynomials, evaluation/substitution, differentiation and division.
 
-* [TypedPolynomials](https://github.com/rdeits/TypedPolynomials.jl)
-* [DynamicPolynomials](https://github.com/blegat/DynamicPolynomials.jl)
+The following packages provide representations of 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
+
+The following packages extend the interface and/or implement algorithms using the interface:
+
+* [SemialgebraicSets](https://github.com/blegat/SemialgebraicSets.jl) : Sets defined by inequalities and equalities between polynomials and algorithms for solving polynomial systems of equations.
+* [HomotopyContinuation](https://github.com/saschatimme/HomotopyContinuation.jl) : Solving systems of polynomials via homotopy continuation.
+* [MultivariateMoments](https://github.com/blegat/MultivariateMoments.jl) : Moments of multivariate measures and their scalar product with polynomials.
+* [PolyJuMP](https://github.com/JuliaOpt/PolyJuMP.jl) : A [JuMP](https://github.com/JuliaOpt/JuMP.jl) extension for Polynomial Optimization.
+* [SumOfSquares](https://github.com/JuliaOpt/SumOfSquares.jl) : Certifying the nonnegativity of polynomials, minimizing/maximizing polynomials and optimization over sum of squares polynomials using Sum of Squares Programming.
+
+## Contents
 ```@contents
-Pages = ["apireference.md"]
+Pages = ["types.md", "substitution.md", "differentiation.md", "division.md"]
 Depth = 3
 ```
diff --git a/docs/src/substitution.md b/docs/src/substitution.md
@@ -0,0 +1,15 @@
+# Subtitution
+
+Given a polynomial, say ``p(x, y) = 3x^2y + x + 2y + 1``, one can evaluate it at a given point, e.g. ``p(2, 1) = 12 + 2 + 2 + 1 = 17`` or substitute one or more variable by a value or polynomial, e.g. ``p(x, xy^2 + 1) = 3x^2(xy^2+1) + x + 2(xy^2+1) + 1 = 3x^3y^2 + 2xy^2 + 3x^2 + x + 3``.
+We distinguish the two operation as follows
+
+* We call an evaluation an operation where **every** variable should be replace by a new value or polynomial, the syntax is `p(x => 2, y => 1)`.
+* We call a subsitution an operation where **some** (or all variables) are subtituted into a new value or polynomial, the syntax is `subs(p, y => x*y^2 + 1)`.
+
+The distinction is important for type stability for some implementations (it is important for [DynamicPolynomials](https://github.com/blegat/DynamicPolynomials.jl) but not for [TypedPolynomials](https://github.com/rdeits/TypedPolynomials.jl)).
+Indeed consider a polynomial with `Int` coefficients for which we ask to replace some variables with `Int` values. If all the variables are replaced with `Int`s, the return type should be `Int`.
+However, if some variables only are replaced by `Int` then the return type should be a polynomial with `Int` coefficients.
+
+```@docs
+subs
+```
diff --git a/docs/src/apireference.md → docs/src/types.md b/docs/src/apireference.md → docs/src/types.md
@@ -2,11 +2,12 @@
 CurrentModule = MultivariatePolynomials
 ```
 
-# API
+# Types
 
 ## Variables
 
 ```@docs
+AbstractVariable
 variable
 name
 similarvariable
@@ -16,21 +17,24 @@ similarvariable
 ## Monomials
 
 ```@docs
+AbstractMonomialLike
+AbstractMonomial
 monomialtype
 variables
 nvariables
 exponents
 degree
 isconstant
 powers
-divides
 constantmonomial
 mapexponents
 ```
 
 ## Terms
 
 ```@docs
+AbstractTermLike
+AbstractTerm
 term
 termtype
 coefficient
@@ -43,6 +47,8 @@ zeroterm
 ## Polynomials
 
 ```@docs
+AbstractPolynomialLike
+AbstractPolynomial
 polynomial
 polynomialtype
 terms
@@ -60,6 +66,11 @@ removemonomials
 monic
 ```
 
+## Rational Polynomial Function
+
+A rational polynomial function can be constructed with the `/` operator. Common operations such as `+`, `-`, `*`, `-` have been implemented between rational functions.
+The numerator and denominator polynomials can be retrieved by the `numerator` and `denominator` functions.
+
 ## Monomial Vectors
 
 ```@docs

diff --git a/src/MultivariatePolynomials.jl b/src/MultivariatePolynomials.jl
@@ -9,38 +9,77 @@ import Base: *, +, -, /, ^, ==,
     one, zero, transpose, isapprox, @pure, dot, copy
 
 export AbstractPolynomialLike, AbstractTermLike, AbstractMonomialLike
+"""
+    AbstractPolynomialLike{T}
+
+Abstract type for a value that can act like a polynomial. For instance, an `AbstractTerm{T}` is an `AbstractPolynomialType{T}` since it can act as a polynomial of only one term.
+"""
 abstract type AbstractPolynomialLike{T} end
+"""
+    AbstractTermLike{T}
+
+Abstract type for a value that can act like a term. For instance, an `AbstractMonomial` is an `AbstractTermLike{Int}` since it can act as a term with coefficient `1`.
+"""
 abstract type AbstractTermLike{T} <: AbstractPolynomialLike{T} end
+"""
+    AbstractMonomialLike
+
+Abstract type for a value that can act like a monomial. For instance, an `AbstractVariable` is an `AbstractMonomialLike` since it can act as a monomial of one variable with degree `1`.
+"""
 abstract type AbstractMonomialLike <: AbstractTermLike{Int} end
 
 export AbstractVariable, AbstractMonomial, AbstractTerm, AbstractPolynomial
+"""
+    AbstractVariable <: AbstractMonomialLike
+
+Abstract type for a variable.
+"""
 abstract type AbstractVariable <: AbstractMonomialLike end
+
+"""
+    AbstractMonomial <: AbstractMonomialLike
+
+Abstract type for a monomial, i.e. a product of variables elevated to a nonnegative integer power.
+"""
 abstract type AbstractMonomial <: AbstractMonomialLike end
+
+"""
+    AbstractTerm{T} <: AbstractTerm{T}
+
+Abstract type for a term of coefficient type `T`, i.e. the product between a value of type `T` and a monomial.
+"""
 abstract type AbstractTerm{T} <: AbstractTermLike{T} end
+
+"""
+    AbstractPolynomial{T} <: AbstractPolynomialLike{T}
+
+Abstract type for a polynomial of coefficient type `T`, i.e. a sum of `AbstractTerm{T}`s.
+"""
 abstract type AbstractPolynomial{T} <: AbstractPolynomialLike{T} end
 
 const APL{T} = AbstractPolynomialLike{T}
 
 include("zip.jl")
+
 include("variable.jl")
 include("monomial.jl")
 include("term.jl")
-include("monovec.jl")
 include("polynomial.jl")
+include("monovec.jl")
 
 include("rational.jl")
 
-include("conversion.jl")
-include("promote.jl")
+include("show.jl")
 include("hash.jl")
 
+include("promote.jl")
+include("conversion.jl")
+
 include("operators.jl")
-include("comparision.jl")
+include("comparison.jl")
 
 include("substitution.jl")
 include("differentiation.jl")
 include("division.jl")
 
-include("show.jl")
-
 end # module
diff --git a/src/comparision.jl → src/comparison.jl b/src/comparision.jl → src/comparison.jl
diff --git a/src/differentiation.jl b/src/differentiation.jl
@@ -1,6 +1,29 @@
 # I do not use it but I import the function to add a method
 export differentiate
 
+"""
+    differentiate(p::AbstractPolynomialLike, v::AbstractVariable, deg::Int=1)
+
+Differentiate `deg` times the polynomial `p` by the variable `v`.
+
+    differentiate(p::AbstractPolynomialLike, vs, deg::Int=1)
+
+Differentiate `deg` times the polynomial `p` by the variables of the vector or tuple of variable `vs` and return an array of dimension `deg`.
+
+    differentiate(p::AbstractArray{<:AbstractPolynomialLike, N}, vs, deg::Int=1) where N
+
+Differentiate the polynomials in `p` by the variables of the vector or tuple of variable `vs` and return an array of dimension `N+deg`.
+
+### Examples
+
+```julia
+p = 3x^2*y + x + 2y + 1
+differentiate(p, x) # should return 6xy + 1
+differentiate(p, (x, y)) # should return [6xy+1, 3x^2+1]
+```
+"""
+function differentiate end
+
 # Fallback for everything else
 _diff_promote_op(::Type{T}, ::Type{<:AbstractVariable}) where T = T
 differentiate(α::T, v::AbstractVariable) where T = zero(T)

diff --git a/src/division.jl b/src/division.jl
@@ -1,13 +1,14 @@
-export monic
+export divides
 
 """
     divides(t1::AbstractTermLike, t2::AbstractTermLike)
 
 Returns whether the monomial of t1 divides the monomial of t2.
 
-# Examples
+### Examples
 
-Calling `divides(3x*y, 2x^2*y)` should return true but calling `divides(x*y, y)` should return false.
+Calling `divides(2x^2y, 3xy)` should return false because `x^2y` does not divide `xy` since `x` has a degree 2 in `x^2y` which is greater than the degree of `x` on `xy`.
+However, calling `divides(3xy, 2x^2y)` should return true.
 """
 function divides(t1::AbstractTermLike, t2::AbstractTermLike)
     divides(monomial(t1), monomial(t2))

diff --git a/src/monomial.jl b/src/monomial.jl
@@ -1,4 +1,4 @@
-export variables, nvariables, exponents, degree, isconstant, powers, divides, constantmonomial, mapexponents
+export variables, nvariables, exponents, degree, isconstant, powers, constantmonomial, mapexponents
 
 """
     monomialtype(p::AbstractPolynomialLike)
@@ -11,7 +11,7 @@ Returns the type of the monomials of a polynomial of type `PT`.
 """
 monomialtype(::Union{M, Type{M}}) where M<:AbstractMonomial = M
 monomialtype(::Union{PT, Type{PT}}) where PT <: APL = monomialtype(termtype(PT))
-monomialtype(::Union{AbstractVector{PT}, Type{AbstractVector{PT}}}) where PT <: APL = monomialtype(PT)
+monomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = monomialtype(PT)
 
 """
     variables(p::AbstractPolynomialLike)
@@ -99,18 +99,6 @@ Calling `powers(3x^4*y) should return `((x, 4), (y, 1))`.
 """
 powers(t::AbstractTermLike) = tuplezip(variables(t), exponents(t))
 
-"""
-    divides(t1::AbstractTermLike, t2::AbstractTermLike)
-
-Returns whether `monomial(t1)` divides `monomial(t2)`.
-
-### Examples
-
-Calling `divides(2x^2y, 3xy)` should return false because `x^2y` does not divide `xy` since `x` has a degree 2 in `x^2y` which is greater than the degree of `x` on `xy`.
-However, calling `divides(3xy, 2x^2y)` should return true.
-"""
-function divides end
-
 """
     constantmonomial(p::AbstractPolynomialType)
 

diff --git a/src/polynomial.jl b/src/polynomial.jl
@@ -2,7 +2,7 @@ export polynomial, polynomialtype, terms, nterms, coefficients, monomials
 export coefficienttype, monomialtype
 export mindegree, maxdegree, extdegree
 export leadingterm, leadingcoefficient, leadingmonomial
-export removeleadingterm, removemonomials
+export removeleadingterm, removemonomials, monic
 
 Base.norm(p::AbstractPolynomialLike, r::Int=2) = norm(coefficients(p), r)
 
@@ -60,11 +60,30 @@ function polynomial(Q::AbstractMatrix, mv::AbstractVector, ::Type{T}) where T
     polynomial(polynomial(Q, mv), T)
 end
 
+"""
+    polynomialtype(p::AbstractPolynomialLike)
+
+Returns the type that `p` would have if it was converted into a polynomial.
+
+    polynomialtype(::Type{PT}) where PT<:AbstractPolynomialLike
+
+Returns the same as `polynomialtype(::PT)`.
+
+    polynomialtype(p::AbstractPolynomialLike, ::Type{T}) where T
+
+Returns the type that `p` would have if it was converted into a polynomial of coefficient type `T`.
+
+    termtype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomialLike, T}
+
+Returns the same as `polynomialtype(::PT, ::Type{T})`.
+"""
 polynomialtype(::Union{P, Type{P}}) where P <: APL = Base.promote_op(polynomial, P)
 polynomialtype(::Union{P, Type{P}}) where P <: AbstractPolynomial = P
 polynomialtype(::Union{M, Type{M}}) where M<:AbstractMonomialLike = polynomialtype(termtype(M))
 polynomialtype(::Union{M, Type{M}}, ::Type{T}) where {M<:AbstractMonomialLike, T} = polynomialtype(termtype(M, T))
 polynomialtype(::Union{P, Type{P}}, ::Type{T}) where {P <: APL, T} = polynomialtype(polynomialtype(P), T)
+polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = polynomialtype(PT)
+polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}, ::Type{T}) where {PT <: APL, T} = polynomialtype(PT, T)
 
 function uniqterms(ts::AbstractVector{T}) where T <: AbstractTerm
     result = T[]
@@ -308,6 +327,8 @@ end
 
 """
     monic(p::AbstractPolynomialLike)
+
+Returns `p / leadingcoefficient(p)` where the leading coefficient of the returned polynomials is made sure to be exactly one to avoid rounding error.
 """
 function monic(p::APL)
     α = leadingcoefficient(p)

diff --git a/src/rational.jl b/src/rational.jl
@@ -7,6 +7,9 @@ struct RationalPoly{NT <: APL, DT <: APL}
     den::DT
 end
 
+Base.numerator(r::RationalPoly) = r.num
+Base.denominator(r::RationalPoly) = r.den
+
 Base.convert(::Type{RationalPoly{NT, DT}}, q::RationalPoly{NT, DT}) where {NT, DT} = q
 Base.convert(::Type{RationalPoly{NTout, DTout}}, q::RationalPoly{NTin, DTin}) where {NTout, DTout, NTin, DTin} = convert(NTout, q.num) / convert(DTout, q.den)
 #function Base.convert(::Type{RationalPoly{NT, DT}}, p::NT) where {NT, DT}