Merge pull request #274 from rdbliss/resultant

Add resultant methods to fmpz_poly and fmpq_poly

Author: oscarbenjamin

README.md

@@ -162,6 +162,18 @@ CHANGELOG
 Next release (0.8.0)...
 -----------------------
 
+Contributors
+
+- Robert Dougherty-Bliss (RDB)
+
+Changes
+
+- [gh-274](https://github.com/flintlib/python-flint/pull/274),
+  Add resultant methods to `fmpz_poly`, `fmpq_poly` and
+  `nmod_poly`. Now all univariate and polynomial types have the
+  resultant method except for `fq_default_poly`. (RDB)
+
+
 0.7.0
 -----
 

src/flint/test/test_all.py

@@ -2834,6 +2834,42 @@ def setbad(obj, i, val):
         if type(p) == flint.fq_default_poly:
             assert raises(lambda: p.integral(), NotImplementedError)
 
+        # resultant checks.
+        x = P([0, 1])
+
+        if composite_characteristic and type(x) in [flint.fmpz_mod_poly, flint.nmod_poly]:
+            # Flint sometimes crashes in this case, even though the resultant
+            # could be computed.
+            divisor = characteristic.factor()[0][0]
+            assert raises(lambda: x.resultant(x + divisor), ValueError)
+        elif type(x) == flint.fq_default_poly:
+            # Flint does not implement resultants over GF(q) for nonprime q, so
+            # there's nothing for us to check.
+            pass
+        else:
+            assert x.resultant(x) == 0
+            assert x.resultant(x**2 + x - x) == 0
+            assert x.resultant(x**10 - x**5 + 1) == S(1)
+            assert (x - 1).resultant(x**5 + 1) == S(2)
+
+            for k in range(-10, 10):
+                assert x.resultant(x + S(k)) == S(k)
+
+def test_poly_resultants():
+    # Check that the resultant of two cyclotomic polynomials is right.
+    # See Dresden's 2012 "Resultants of Cyclotomic Polynomials"
+    for m in range(1, 50):
+        for n in range(m + 1, 50):
+            a = flint.fmpz_poly.cyclotomic(m)
+            b = flint.fmpz_poly.cyclotomic(n)
+            q, r = divmod(flint.fmpz(n), flint.fmpz(m))
+            fs = q.factor()
+            if r != 0 or len(fs) > 1:
+                assert a.resultant(b) == 1
+            else:
+                prime = fs[0][0]
+                tot = flint.fmpz(m).euler_phi()
+                assert a.resultant(b) == prime**tot
 
 def _all_mpolys():
     return [
@@ -2869,7 +2905,6 @@ def _all_mpolys():
         ),
     ]
 
-
 def test_mpolys():
     for P, get_context, S, is_field, characteristic in _all_mpolys():
 
@@ -4832,6 +4867,8 @@ def test_all_tests():
     test_polys,
     test_mpolys,
 
+    test_poly_resultants,
+
     test_fmpz_mpoly_vec,
 
     test_matrices_eq,

src/flint/types/fmpq_poly.pyx

@@ -415,6 +415,32 @@ cdef class fmpq_poly(flint_poly):
         fmpq_poly_gcd(res.val, self.val, (<fmpq_poly>other).val)
         return res
 
+    def resultant(self, other):
+        """
+        Returns the resultant of *self* and *other*.
+
+            >>> A = fmpq_poly([1, 0, -1]); B = fmpq_poly([1, -1])
+            >>> A.resultant(B)
+            0
+            >>> C = fmpq_poly([1, 0, 0, 0, 0, -1, 1])
+            >>> D = fmpq_poly([1, 0, 0, -1, 0, 0, 1])
+            >>> C.resultant(D)
+            3
+            >>> f = fmpq_poly([1, -1] + [0] * 98 + [1])
+            >>> g = fmpq_poly([1] + [0] * 50 + [-1] + [0] * 48 + [1])
+            >>> f.resultant(g)
+            1125899906842623
+
+        """
+        cdef fmpq res
+        other = any_as_fmpq_poly(other)
+        if other is NotImplemented:
+            raise TypeError("cannot convert input to fmpq_poly")
+
+        res = fmpq.__new__(fmpq)
+        fmpq_poly_resultant(res.val, self.val, (<fmpq_poly>other).val)
+        return res
+
     def xgcd(self, other):
         cdef fmpq_poly res1, res2, res3
         other = any_as_fmpq_poly(other)

src/flint/types/fmpz_mod_poly.pyx

@@ -1465,6 +1465,9 @@ cdef class fmpz_mod_poly(flint_poly):
         """
         cdef fmpz_mod res
 
+        if not self.ctx.mod.is_prime():
+            raise ValueError("cannot compute fmpz_mod_poly resultants with composite moduli")
+
         other = self.ctx.any_as_fmpz_mod_poly(other)
         if other is NotImplemented:
             raise TypeError(f"Cannot interpret {other} as a polynomial")

src/flint/types/fmpz_poly.pyx

@@ -397,6 +397,32 @@ cdef class fmpz_poly(flint_poly):
         fmpz_poly_gcd(res.val, self.val, (<fmpz_poly>other).val)
         return res
 
+    def resultant(self, other):
+        """
+        Returns the resultant of *self* and *other*.
+
+            >>> A = fmpz_poly([1, 0, -1]); B = fmpz_poly([1, -1])
+            >>> A.resultant(B)
+            0
+            >>> C = fmpz_poly([1, 0, 0, 0, 0, -1, 1])
+            >>> D = fmpz_poly([1, 0, 0, -1, 0, 0, 1])
+            >>> C.resultant(D)
+            3
+            >>> f = fmpz_poly([1, -1] + [0] * 98 + [1])
+            >>> g = fmpz_poly([1] + [0] * 50 + [-1] + [0] * 48 + [1])
+            >>> f.resultant(g)
+            1125899906842623
+
+        """
+        cdef fmpz res
+        other = any_as_fmpz_poly(other)
+        if other is NotImplemented:
+            raise TypeError("cannot convert input to fmpz_poly")
+
+        res = fmpz.__new__(fmpz)
+        fmpz_poly_resultant(res.val, self.val, (<fmpz_poly>other).val)
+        return res
+
     def factor(self):
         """
         Factors self into irreducible factors, returning a tuple

src/flint/types/nmod_poly.pyx

@@ -621,6 +621,31 @@ cdef class nmod_poly(flint_poly):
         nmod_poly_gcd(res.val, self.val, (<nmod_poly>other).val)
         return res
 
+    def resultant(self, other):
+        """
+        Returns the resultant of *self* and *other*.
+
+            >>> f = nmod_poly([1, 2, 3], 3)
+            >>> g = nmod_poly([1, 0, 1], 3)
+            >>> f.resultant(f)
+            0
+            >>> f.resultant(g)
+            2
+
+        """
+        cdef ulong res
+
+        mod = any_as_fmpz(self.val.mod.n)
+        if not mod.is_prime():
+            raise ValueError("cannot compute nmod_poly resultants with composite moduli")
+
+        other = any_as_nmod_poly(other, (<nmod_poly>self).val.mod)
+        if other is NotImplemented:
+            raise TypeError("cannot convert input to nmod_poly")
+
+        res = nmod_poly_resultant(self.val, (<nmod_poly>other).val)
+        return res
+
     def xgcd(self, other):
         cdef nmod_poly res1, res2, res3
         other = any_as_nmod_poly(other, (<nmod_poly>self).val.mod)