feat: add resultant for list of multivariate polynomials

src/Rings/Rings.jl

@@ -36,3 +36,4 @@ include("PBWAlgebraQuo.jl")
 include("FreeAssAlgIdeal.jl")
 include("hilbert.jl")
 include("primary_decomposition_helpers.jl")
+include("resultant.jl")

src/Rings/resultant.jl

@@ -0,0 +1,156 @@
+################################################################################
+#
+#  Classical resultant
+#
+################################################################################
+
+# Classical resultant for families of polynomials. Inspired by
+# "A package for computations with classical resultants"
+# by Stagliano, http://dx.doi.org/10.2140/jsag.2018.8.21
+
+@doc raw"""
+    resultant(F::Vector{MPolyRingElem})
+
+Return the Macaulay resultant of a list of $n$ homogeneous polynomials in $n$
+variables.
+
+# Examples
+
+```jldoctest
+julia> QQt, t = QQ["t"];
+
+julia> R, (x1, x2, x0) = QQt["x1", "x2", "x0"];
+
+julia> F = [x2^2 - t*x0*x1 , x2^2 - t*x0*x2, x1^2 + x2^2 - x0^2];
+
+julia> resultant(F)
+-2*t^6 + t^4
+```
+"""
+function resultant(F::Vector{<: MPolyRingElem})
+  # top-level resultant function, which does the argument verification
+  @req length(F) > 0 "Resultant not defined for zero polynomials"
+  R = parent(F[1])
+  n = ngens(R) - 1
+  @req n + 1 == length(F) "Number of polynomials ($(length(F))) must be number of variables ($(n + 1))"
+  # Cannot test the following, because it is not implemented for all
+  # rings
+  # @req all(is_homogeneous, F) "Polynomials must be homogeneous)
+  return _resultant(F)
+end
+
+function _resultant(F::Vector{<: MPolyRingElem{<: FieldElem}})
+  # field case, apply Poisson formula directly
+  return _resultant_poisson(F)
+end
+
+function _resultant(F::Vector{<: MPolyRingElem})
+  # ring case, pass to the field of fractions if possible
+  R = coefficient_ring(parent(F[1]))
+  if is_domain_type(elem_type(R))
+    K = fraction_field(R)
+    Kx, x = polynomial_ring(K, ngens(parent(F[1])); cached = false)
+    FK = map_coefficients.(Ref(K), F, cached = false, parent = Kx)
+    r = _resultant_poisson(FK)
+    @assert is_one(denominator(r))
+    return numerator(r)
+  else
+    error("Resultant implemented only for domains. Please open an issue")
+  end
+end
+
+function _resultant_poisson(F::Vector{<: MPolyRingElem{<:FieldElem}})
+  # Resultant via Poisson formula apres
+  # Theorem 3.4, p. 96 of Cox-Little-O'shea, Using Algebraic Geometry (2005)
+  #
+  # In case a naive application returns zero, it is inconclusive.
+  # In this case we compute the dimension, followed by a random coordinate
+  # transformation in the zero-dimensional case.
+
+  R = parent(F[1])
+  r = __resultant_poisson(F)
+  if !is_zero(r)
+    return r
+  end
+
+  # now r == 0, but we could have had bad projections
+  d = dim(ideal(F))
+
+  if d > 0
+    # there are infinitely many solutions, hence a nontrivial one
+    return r
+  end
+
+  # we make a random SL_n coordinate change, which does not change the resultant
+  # See 1.5 (and the references) of "A package for computations with classical
+  # resultants" by Stagliano, http://dx.doi.org/10.2140/jsag.2018.8.21
+
+  # keep a counter to capture bugs
+  k = 0
+  while is_zero(r)
+    k > 100 && error("Something wrong in the resultant. Please report this bug")
+    h = _random_sln_transformation(R)
+    r = __resultant_poisson(h.(F))
+  end
+  return r
+end
+
+function __resultant_poisson(F::Vector{<: MPolyRingElem{<:FieldElem}})
+  # If we hit a zero, it is inconclusive and we return the zero
+  #@req all(is_homogeneous, F) "Polynomials must be homogenous"
+  R = parent(F[1])
+  n = ngens(R) - 1
+  K = coefficient_ring(R)
+  d = total_degree.(F)
+  # Could have more special cases
+  if n == 0 # one polynomial
+    return is_zero(F[1]) ? zero(K) : leading_coefficient(F[1])
+  end
+  S, x = polynomial_ring(K, :x => 0:n-1, cached = false)
+  h = hom(R, S, vcat(x, [one(S)]))
+  f = h.(F)
+  h = hom(R, S, vcat(x, [zero(S)]))
+  Fbar = h.(F)
+  res0 = __resultant_poisson(Fbar[1:end-1])
+  if is_zero(res0)
+    return res0
+  end
+  I = ideal(S, f[1:end-1])
+  if dim(I) > 0
+    return zero(K)
+  end
+  Q, mQ = quo(S, I)
+  V, VtoQ = vector_space(K, Q)
+  @assert dim(V) == prod(d[1:end-1]; init = 1)
+  M = zero_matrix(K, dim(V), dim(V))
+  fn = mQ(last(f))
+  for i in 1:dim(V)
+    v = VtoQ\(VtoQ(V[i]) * fn)
+    for j in 1:dim(V)
+      M[i, j] = v[j]
+    end
+  end
+  return res0^d[end] * det(M)
+end
+
+function _random_sln_transformation(R::MPolyRing)
+  K = base_ring(R)
+  n = ngens(R)
+  @assert K isa Field
+  A = identity_matrix(K, n)
+  if n <= 1
+    return A
+  end
+  # do five random elementary operations
+  for l in 1:5
+    i = rand(1:n)
+    j = rand(1:n)
+    while j == i
+      j = rand(1:n)
+    end
+    add_row!(A, one(K), i, j)
+  end
+  X = gens(R)
+  h = hom(R, R, [sum((A[i, j] * X[j] for j in 1:n); init = zero(R)) for i in 1:n], check = false)
+  return h
+end

test/Rings/resultant.jl

@@ -0,0 +1,27 @@
+@testset "resultant" begin
+  QQt, t = QQ["t"]
+  R, (x1, x2, x0) = QQt["x1", "x2", "x0"]
+  F = [x2^2 - t*x0*x1 , x2^2-t*x0*x2, x1^2+x2^2-x0^2]
+  # this example requires coordinate change
+  @test resultant(F) == -2*t^6 + t^4
+
+  R, (x, y, z) = QQ["x", "y", "z"]
+  F0 = x^3 + y^2 * z
+  F1 = x*y + y^2 + x*z + y * z
+  F2 = y^4 + z^4
+  r = resultant([F0, F1, F2])
+  @test parent(r) === QQ && r == 16
+
+  @test_throws ArgumentError resultant(typeof(F0)[])
+  @test_throws ArgumentError resultant([F0])
+  @test_throws ArgumentError resultant([F0, F0, F0, F0])
+
+  # example from the Macaulay2 documentation:
+   QQab, (a, b) = QQ[:a, :b]
+   R, (x, y, z, w) = QQab[:x, :y, :z, :w]
+   F = [(7//3)*x+(7//2)*y+z+2*w,
+        ((10//7)*a+b)*x^2+(a+(5//4)*b)*x*y+(2*a+(1//2)*b)*y^2+((7//8)*a+(7//5)*b)*x*z+((3//4)*a+b)*y*z+((7//8)*a+(1//7)*b)*z^2+((5//7)*a+(4//3)*b)*x*w+(9*a+10*b)*y*w+((7//5)*a+(3//4)*b)*z*w+((4//3)*a+5*b)*w^2,
+        ((1//2)*a+(7//5)*b)*x^3+((1//2)*a+10*b)*x^2*y+((8//9)*a+(3//5)*b)*x*y^2+(a+(7//6)*b)*y^3+((3//7)*a+(3//4)*b)*x^2*z+((1//3)*a+(9//10)*b)*x*y*z+((9//4)*a+b)*y^2*z+((1//6)*a+(1//5)*b)*x*z^2+(3*a+(5//2)*b)*y*z^2+((5//3)*a+(3//7)*b)*z^3+(a+b)*x^2*w+((4//5)*a+(5//4)*b)*x*y*w+((5//3)*a+(5//8)*b)*y^2*w+((3//2)*a+(1//6)*b)*x*z*w+((1//3)*a+(4//5)*b)*y*z*w+(9*a+(1//3)*b)*z^2*w+((7//3)*a+(5//4)*b)*x*w^2+(a+(3//4)*b)*y*w^2+((9//8)*a+(7//8)*b)*z*w^2+((9//7)*a+2*b)*w^3,
+        2*x+(1//4)*y+(8//3)*z+(4//5)*w]
+  @test resultant(F) == -(21002161660529014459938925799//2222549728809984000000)*a^5-(2085933800619238998825958079203//12700284164628480000000)*a^4*b-(348237304382147063838108483692249//889019891523993600000000)*a^3*b^2-(38379949248928909714532254698073//35278567123968000000000)*a^2*b^3-(1146977327343523453866040839029//1119954511872000000000)*a*b^4-(194441910898734675845094443//895963609497600000)*b^5
+end