Trac #16014: Improvements to discriminant computation

Currently (i.e. on sage 6.1.1), the discriminant computation fails with
a segmenttation fault when the coefficients are too complicated. For
example:

{{{
#!python
PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
PRmu.<mu> = PR[]
E1 = diagonal_matrix(PR, [1, b^2, -b^2])
RotScale = matrix(PR, [[1, -t1, 0], [t1, 1, 0], [0, 0, 1]])
E2 = diagonal_matrix(PR, [1, 1,
1+t1^2])*RotScale.transpose()*E1*RotScale
Transl = matrix(PR, [[1, 0, -x1], [0, 1, -y1], [0, 0, 1]])
E3 = Transl.transpose()*E2*Transl
E4 = E3.subs(t1=t2, x1=x2, y1=y2)
det(mu*E3 + E4).discriminant()
}}}

Apparently the Python interpreter reaches some limit on the depth of the
call stack, since it seems to call `symtable_visit_expr` from
`Python-2.7.6/Python/symtable.c:1191` over and over and over again. But
this ticket here is not about fixing Python. (If you want a ticket for
that, feel free to create one and provide a reference here.)

Instead, one should observe that the above is only a discriminant of a
third degree polynomial. As such, it has an easy and well-known formula
for the discriminant. Plugging the complicated coefficients into this
easy formula leads to the solution rather quickly (considering the
immense size of said solution). So perhaps we should employ this
approach in the `discriminant` method.

I've [https://groups.google.com/d/topic/sage-
support/fABfhHyioa0/discussion discussed this question on sage-support],
and will probably continue the discussion. Right now, I'm unsure about
when to choose a new method over the existing one. Should it be based on
degree, using hardcoded versions of the well-known discriminants for
lower degrees? Or should it rather be based on base ring, thus employing
a computation over some simple generators instead of the full
coefficients if the base ring is a polynomial ring or some other fairly
complex object. Here is a code snippet to illustrate the latter idea:

{{{
#!python
def generic_discriminant(p):
    """
    Compute discriminant of univariate (sparse or dense) polynomial.

    The discriminant is 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 like complicated polynomials. 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.
    """
    pr = p.parent()
    if not pr.ngens() == 1:
        raise ValueError("Must be univariate")
    ex = p.exponents()
    pr1 = PolynomialRing(ZZ, "x", len(ex))
    pr2.<y> = pr1[]
    py = sum(v*y^e for v, e in zip(pr1.gens(), ex))
    d = py.discriminant()
    return d(p.coefficients())
}}}

Personally I think I prefer a case distinction based on base ring, or a
combination. If we do the case distinction based on base ring, then I'd
also value some input as to what rings should make use of this new
approach, and how to check for them. For example, I just noticed that
#15061 discusses issues with discriminants of polynomials over power
series. So I guess those might benefit from this approach as well.

URL: http://trac.sagemath.org/16014
Reported by: gagern
Ticket author(s): Martin von Gagern
Reviewer(s): Peter Bruin

src/sage/rings/polynomial/polynomial_element.pyx

@@ -91,7 +91,10 @@ 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.misc.cachefunc import cached_function
 
 from sage.categories.map cimport Map
 from sage.categories.morphism cimport Morphism
@@ -4891,8 +4894,8 @@ cdef class Polynomial(CommutativeAlgebraElement):
 
         .. note::
 
-           Uses the identity `R_n(f) := (-1)^(n (n-1)/2) R(f,
-           f') a_n^(n-k-2)`, where `n` is the degree of self,
+           Uses the identity `R_n(f) := (-1)^{n (n-1)/2} R(f,
+           f') a_n^{n-k-2}`, where `n` is the degree of self,
            `a_n` is the leading coefficient of self, `f'`
            is the derivative of `f`, and `k` is the degree
            of `f'`. Calls :meth:`.resultant`.
@@ -4971,10 +4974,41 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: R.<s> = PolynomialRing(Qp(2))
             sage: (s^2).discriminant()
             0
+
+        TESTS:
+
+        This was fixed by :trac:`16014`::
+
+            sage: PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
+            sage: PRmu.<mu> = PR[]
+            sage: E1 = diagonal_matrix(PR, [1, b^2, -b^2])
+            sage: M = matrix(PR, [[1,-t1,x1-t1*y1],[t1,1,y1+t1*x1],[0,0,1]])
+            sage: E1 = M.transpose()*E1*M
+            sage: E2 = E1.subs(t1=t2, x1=x2, y1=y2)
+            sage: det(mu*E1 + E2).discriminant().degrees()
+            (24, 12, 12, 8, 8, 8, 8)
+
+        This addresses an issue raised by :trac:`15061`::
+
+            sage: R.<T> = PowerSeriesRing(QQ)
+            sage: F = R([1,1],2)
+            sage: RP.<x> = PolynomialRing(R)
+            sage: P = x^2 - F
+            sage: P.discriminant()
+            4 + 4*T + O(T^2)
         """
+        # Late import to avoid cyclic dependencies:
+        from sage.rings.power_series_ring import is_PowerSeriesRing
         if self.is_zero():
             return self.parent().zero_element()
         n = self.degree()
+        base_ring = self.parent().base_ring()
+        if (is_MPolynomialRing(base_ring) or
+            is_PowerSeriesRing(base_ring)):
+            # It is often cheaper to compute discriminant of simple
+            # multivariate polynomial and substitute the real
+            # coefficients into that result (see #16014).
+            return universal_discriminant(n)(list(self))
         d = self.derivative()
         k = d.degree()
 
@@ -7354,6 +7388,42 @@ cdef class Polynomial_generic_dense(Polynomial):
 def make_generic_polynomial(parent, coeffs):
     return parent(coeffs)
 
+@cached_function
+def universal_discriminant(n):
+    r"""
+    Return the discriminant of the 'universal' univariate polynomial
+    `a_n x^n + \cdots + a_1 x + a_0` in `\ZZ[a_0, \ldots, a_n][x]`.
+
+    INPUT:
+
+    - ``n`` - degree of the polynomial
+
+    OUTPUT:
+
+    The discriminant as a polynomial in `n + 1` variables over `\ZZ`.
+    The result will be cached, so subsequent computations of
+    discriminants of the same degree will be faster.
+
+    EXAMPLES::
+
+        sage: from sage.rings.polynomial.polynomial_element import universal_discriminant
+        sage: universal_discriminant(1)
+        1
+        sage: universal_discriminant(2)
+        a1^2 - 4*a0*a2
+        sage: universal_discriminant(3)
+        a1^2*a2^2 - 4*a0*a2^3 - 4*a1^3*a3 + 18*a0*a1*a2*a3 - 27*a0^2*a3^2
+        sage: universal_discriminant(4).degrees()
+        (3, 4, 4, 4, 3)
+
+    .. SEEALSO::
+        :meth:`Polynomial.discriminant`
+    """
+    pr1 = PolynomialRing(ZZ, n + 1, 'a')
+    pr2 = PolynomialRing(pr1, 'x')
+    p = pr2(list(pr1.gens()))
+    return (1 - (n&2))*p.resultant(p.derivative())//pr1.gen(n)
+
 cdef class ConstantPolynomialSection(Map):
     """
     This class is used for conversion from a polynomial ring to its base ring.