Further mpoly interface (#228)

* Add functions for taking (multivariate) coeff of monomial.

* Change terms and monomials functions to work with any MPolyElem.

* Fix typo in comment.

src/AbstractAlgebra.jl

@@ -164,7 +164,8 @@ import .Generic: add!, addeq!, addmul!, base_ring, cached,
                  sort_terms!, sub, subst, swap_rows,
                  swap_rows!, symbols, term, terms, total_degree, to_univariate,
                  trail, truncate, typed_hcat, typed_hvcat, upscale,
-                 valuation, var, vars, weak_popov, weak_popov_with_trafo, zero,
+                 valuation, var, var_index, vars, weak_popov,
+                 weak_popov_with_trafo, zero,
                  zero!, zero_matrix, @PolynomialRing
 
 # Do not export exp, sqrt, numerator and denominator as we define our own
@@ -227,7 +228,7 @@ export add!, addeq!, addmul!, addmul_delayed_reduction!, base_ring, cached,
                  sort_terms!, sub, subst, swap_rows, swap_rows!,
                  term, terms, to_univariate,
                  total_degree, tr, trail, truncate, typed_hcat, typed_hvcat,
-                 upscale, valuation, var, vars,
+                 upscale, valuation, var, var_index, vars,
                  weak_popov, weak_popov_with_trafo, zero, zero!,
                  zero_matrix, @PolynomialRing
 

src/generic/MPoly.jl

@@ -12,7 +12,7 @@ export max_fields, total_degree, gens, divides, isconstant, isdegree,
        derivative, change_base_ring,  to_univariate, degrees, deflation,
        combine_like_terms!, exponent, exponent_vector, exponent_vectors,
        set_exponent_vector!, sort_terms!, coeffs, monomial, monomial!,
-       monomials, term, terms, @PolynomialRing
+       monomials, term, terms, var_index, @PolynomialRing
 
 ###############################################################################
 #
@@ -155,6 +155,29 @@ function vars(p::AbstractAlgebra.Generic.MPoly{T}) where {T <: RingElement}
    return(vars_in_p)
 end
 
+@doc Markdown.doc"""
+   var_index(p::AbstractAlgebra.MPolyElem{T}) where {T <: RingElement}
+> Return the index of the given variable $x$. If $x$ is not a variable
+> in a multivariate polynomial ring, an exception is raised.
+"""
+function var_index(x::AbstractAlgebra.MPolyElem{T}) where {T <: RingElement}
+   !ismonomial(x) && error("Not a variable in var_index")
+   exps = exponent_vector(x, 1)
+   count = 0
+   index = 0
+   for i = 1:length(exps)
+      if exps[i] > 1
+         error("Not a variable in var_index")
+      end
+      if exps[i] == 1
+         count += 1
+         index = i
+      end
+   end
+   count != 1 && error("Not a variable in var_index")
+   return index
+end
+
 @doc Markdown.doc"""
     ordering(a::MPolyRing{T}) where {T <: RingElement}
 > Return the ordering of the given polynomial ring as a symbol. The options are
@@ -169,6 +192,67 @@ function check_parent(a::MPoly{T}, b::MPoly{T}) where T <: RingElement
       error("Incompatible polynomial rings in polynomial operation")
 end
 
+###############################################################################
+#
+#   Multivariate coefficients
+#
+###############################################################################
+
+@doc Markdown.doc"""
+    coeff(a::MPoly{T}, vars::Vector{Int}, exps::Vector{Int}) where T <: RingElement
+> Return the "coefficient" of $a$ (as a multivariate polynomial in the same
+> ring) of the monomial consisting of the product of the variables of the given
+> indices raised to the given exponents (note that not all variables need to
+> appear and the exponents can be zero). E.g. `coeff(f, [1, 3], [0, 2])` returns
+> the coefficient of $x^0*z^2$ in the polynomial $f$ (assuming variables
+> $x, y, z$ in that order).
+"""
+function coeff(a::MPoly{T}, vars::Vector{Int}, exps::Vector{Int}) where T <: RingElement
+   unique(vars) != vars && error("Variables not unique")
+   length(vars) != length(exps) &&
+       error("Number of variables does not match number of exponents")
+   for i = 1:length(vars)
+      if vars[i] < 1 || vars[i] > nvars(parent(a))
+         error("Variable index not in range")
+      end
+      if exps[i] < 0
+         error("Exponent cannot be negative")
+      end
+   end
+   exp_vecs = Array{Vector{Int}, 1}(undef, 0)
+   coeffs = Array{T, 1}(undef, 0)
+   for i = 1:length(a)
+      v = exponent_vector(a, i)
+      flag = true
+      for j = 1:length(vars)
+         if v[vars[j]] != exps[j]
+            flag = false
+            break
+         else
+            v[vars[j]] = 0
+         end 
+      end
+      if flag
+         push!(exp_vecs, v)
+         push!(coeffs, coeff(a, i))
+      end
+   end
+   return parent(a)(coeffs, exp_vecs)
+end
+
+@doc Markdown.doc"""
+    coeff(a::T, vars::Vector{T}, exps::Vector{Int}) where T <: AbstractAlgebra.MPolyElem
+> Return the "coefficient" of $a$ (as a multivariate polynomial in the same
+> ring) of the monomial consisting of the product of the given variables
+> to the given exponents (note that not all variables need to appear and the
+> exponents can be zero). E.g. `coeff(f, [x, z], [0, 2])` returns the
+> coefficient of $x^0*z^2$ in the polynomial $f$. 
+"""
+function coeff(a::T, vars::Vector{T}, exps::Vector{Int}) where T <: AbstractAlgebra.MPolyElem
+   varidx = [var_index(x) for x in vars]
+   return coeff(a, varidx, exps)
+end
+
 ###############################################################################
 #
 #   Monomial operations
@@ -493,11 +577,11 @@ function monomial!(m::MPoly{T}, x::MPoly{T}, i::Int) where T <: RingElement
 end
 
 @doc Markdown.doc"""
-    monomials(x::MPoly)
+    monomials(x::AbstractAlgebebra.MPolyElem)
 > Return an array of the monomials of the nonzero terms of the given
 > polynomial, starting with the most significant term.
 """
-function monomials(x::MPolyElem)
+function monomials(x::AbstractAlgebra.MPolyElem)
    return [monomial(x, i) for i = 1:length(x)]
 end
 
@@ -514,11 +598,11 @@ function term(x::MPoly, i::Int)
 end
 
 @doc Markdown.doc"""
-    terms(x::MPoly)
+    terms(x::AbstractAlgebra.MPolyElem)
 > Return an array of the nonzero terms of the given polynomial, starting with
 > the most significant term.
 """
-function terms(x::MPolyElem)
+function terms(x::AbstractAlgebra.MPolyElem)
    return [term(x, i) for i = 1:length(x)]
 end
 

test/generic/MPoly-test.jl

@@ -80,6 +80,7 @@ function test_gen_mpoly_manipulation()
          @test isgen(g[i])
          @test !isgen(g[i] + 1)
          @test gen(S, i) == g[i]
+         @test var_index(gen(S, i)) == i
       end
 
       nv = rand(1:num_vars)
@@ -186,6 +187,31 @@ function test_gen_mpoly_manipulation()
    println("PASS")
 end
 
+function test_gen_mpoly_multivariate_coeff()
+   print("Generic.MPoly.multivariate_coeff...")
+
+   R = ZZ
+
+   for iter = 1:5
+      ord = rand_ordering()
+
+      S, (x, y, z) = PolynomialRing(R, ["x", "y", "z"]; ordering=ord)
+
+      f = -8*x^5*y^3*z^5+9*x^5*y^2*z^3-8*x^4*y^5*z^4-10*x^4*y^3*z^2+8*x^3*y^2*z-10*x*y^3*
+z^4-4*x*y-10*x*z^2+8*y^2*z^5-9*y^2*z^3
+
+      @test coeff(f, [1], [1]) == -10*y^3*z^4-4*y-10*z^2
+      @test coeff(f, [2, 3], [3, 2]) == -10*x^4
+      @test coeff(f, [1, 3], [4, 5]) == 0
+
+      @test coeff(f, [x], [1]) == -10*y^3*z^4-4*y-10*z^2
+      @test coeff(f, [y, z], [3, 2]) == -10*x^4
+      @test coeff(f, [x, z], [4, 5]) == 0
+   end
+
+   println("PASS")
+end
+
 function test_gen_mpoly_total_degree()
    print("Generic.MPoly.total_degree...")
    max = 50
@@ -883,6 +909,7 @@ end
 function test_gen_mpoly()
    test_gen_mpoly_constructors()
    test_gen_mpoly_manipulation()
+   test_gen_mpoly_multivariate_coeff()
    test_gen_mpoly_unary_ops()
    test_gen_mpoly_binary_ops()
    test_gen_mpoly_adhoc_binary()