From d5e03075e3a9a8beac80396027413256d9b39289 Mon Sep 17 00:00:00 2001
From: Martin von Gagern <Martin.vGagern@gmx.net>
Date: Thu, 27 Mar 2014 00:54:45 +0100
Subject: [PATCH] Use simple coefficients if base is multivariate polynomial.

If the base ring is a multivariate polynomial ring, then the discriminant is
now computed in two steps.  First all non-zero coefficients of the
polynomial are replaced with variables, and the discriminant is computed
using these variables.  Then the original coefficients are substituted into
the result of this computation.  This should be faster in many cases,
particularly with big coefficients.  In those cases, the matrix computation
on simple generators is a lot faster, and the final substitution is cheaper
than carrying all that information along at all times.
---
 .../rings/polynomial/polynomial_element.pyx   | 31 ++++++++++++++-----
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/src/sage/rings/polynomial/polynomial_element.pyx b/src/sage/rings/polynomial/polynomial_element.pyx
index ec433ba0b1d..1624d77b5ee 100644
--- a/src/sage/rings/polynomial/polynomial_element.pyx
+++ b/src/sage/rings/polynomial/polynomial_element.pyx
@@ -88,7 +88,9 @@ import polynomial_fateman
 
 from sage.rings.integer cimport Integer
 from sage.rings.ideal import is_Ideal
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
+from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing
 
 from sage.categories.map cimport Map
 from sage.categories.morphism cimport Morphism
@@ -4990,8 +4992,18 @@ cdef class Polynomial(CommutativeAlgebraElement):
         """
         if self.is_zero():
             return self.parent().zero_element()
-        n = self.degree()
-        d = self.derivative()
+        poly = self
+        if is_MPolynomialRing(self.parent().base_ring()):
+            # It is often cheaper to compute discriminant of simple
+            # multivariate polynomial and substitute the real
+            # coefficients into that result (see #16014).
+            ex = self.exponents()
+            pr1 = PolynomialRing(ZZ, "x", len(ex))
+            pr2 = PolynomialRing(pr1, "y")
+            y = pr2.gen()
+            poly = sum(v*(y**e) for v, e in zip(pr1.gens(), ex))
+        n = poly.degree()
+        d = poly.derivative()
         k = d.degree()
 
         r = n % 4
@@ -4999,17 +5011,20 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if r == 0 or r == 1:
             u = 1
         try:
-            an = self[n]**(n - k - 2)
+            an = poly[n]**(n - k - 2)
         except ZeroDivisionError:
             assert(n-k-2 == -1)
             # Rather than dividing the resultant by the leading coefficient,
             # we alter the Sylvester matrix (see #11782).
-            mat = self.sylvester_matrix(d)
-            mat[0, 0] = self.base_ring()(1)
-            mat[n - 1, 0] = self.base_ring()(n)
-            return u * mat.determinant()
+            mat = poly.sylvester_matrix(d)
+            mat[0, 0] = poly.base_ring()(1)
+            mat[n - 1, 0] = poly.base_ring()(n)
+            res = u * mat.determinant()
         else:
-            return self.base_ring()(u * self.resultant(d) * an)
+            res = poly.base_ring()(u * poly.resultant(d) * an)
+        if poly is not self:
+            res = res(self.coefficients())
+        return res
 
     def reverse(self, degree=None):
         """