Make fmpq_poly.factor() return primitive factors

Author: oscarbenjamin

src/flint/flint_base/flint_base.pyx

@@ -4,6 +4,7 @@ from flint.flintlib.flint cimport (
     __FLINT_RELEASE as _FLINT_RELEASE,
     slong
 )
+from flint.utils.flint_exceptions import DomainError
 from flint.flintlib.mpoly cimport ordering_t
 from flint.flint_base.flint_context cimport thectx
 from flint.flint_base.flint_base cimport Ordering
@@ -249,9 +250,13 @@ cdef class flint_poly(flint_elem):
         roots = []
         factors = self.factor()
         for fac, m in factors[1]:
-            if fac.degree() == fac[1] == 1:
-                v = - fac[0]
-                roots.append((v, m))
+            if fac.degree() == 1:
+                try:
+                    v = - fac[0] / fac[1]
+                except DomainError:
+                    pass
+                else:
+                    roots.append((v, m))
         return roots
 
     def complex_roots(self):

src/flint/test/test_all.py

@@ -3233,6 +3233,85 @@ def test_fmpz_mpoly_vec():
         assert raises(lambda: vec.__setitem__(0, ctx1.from_dict({})), IncompatibleContextError)
 
 
+def _all_polys_mpolys():
+
+    for P, S, is_field in _all_polys():
+        x = P([0, 1])
+        y = None
+        assert isinstance(x, (
+            flint.fmpz_poly,
+            flint.fmpq_poly,
+            flint.nmod_poly,
+            flint.fmpz_mod_poly,
+            flint.fq_default_poly,
+        ))
+        characteristic_zero = isinstance(x, (flint.fmpz_poly, flint.fmpq_poly))
+        yield P, S, [x, y], is_field, characteristic_zero
+
+    for P, ctx_type, S, is_field in _all_mpolys():
+        ctx = ctx_type(2, flint.Ordering.lex, ["x", "y"])
+        x, y = ctx.gens()
+        assert isinstance(x, (
+            flint.fmpz_mpoly,
+            flint.fmpq_mpoly,
+        ))
+        characteristic_zero = isinstance(x, (flint.fmpz_mpoly, flint.fmpq_mpoly))
+        yield P, S, [x, y], is_field, characteristic_zero
+
+
+def test_factor_poly_mpoly():
+    """Test that factor() is consistent across different poly/mpoly types."""
+
+    def factor(p):
+        coeff, factors = p.factor()
+        try:
+            lc = p.leading_coefficient()
+        except AttributeError:
+            # XXX: Not all univariate types have a leading_coefficient method.
+            lc = p[0]
+        assert type(coeff) is type(lc)
+        assert isinstance(factors, list)
+        assert all(isinstance(f, tuple) for f in factors)
+        for fac, m in factors:
+            assert type(fac) is type(p)
+            assert type(m) is int
+        return coeff, sorted(factors, key=lambda p: (p[1], str(p[0])))
+
+    for P, S, [x, y], is_field, characteristic_zero in _all_polys_mpolys():
+
+        assert factor(x) == (S(1), [(x, 1)])
+        assert factor(x**2) == (S(1), [(x, 2)])
+        assert factor(x*(x+1)) == (S(1), [(x, 1), (x+1, 1)])
+        assert factor(2*(x+1)) == (S(2), [(x+1, 1)])
+
+        if characteristic_zero:
+            # primitive factors over Z for Z and Q.
+            assert factor(2*x+1) == (S(1), [(2*x+1, 1)])
+        else:
+            # monic factors over Z/pZ and GF(p^d)
+            assert factor(2*x+1) == (S(2), [(x+S(1)/2, 1)])
+
+        if is_field:
+            if characteristic_zero:
+                # primitive factors over Z for Z and Q.
+                assert factor((2*x+1)/7) == (S(1)/7, [(2*x+1, 1)])
+            else:
+                # monic factors over Z/pZ and GF(p^d)
+                assert factor((2*x+1)/7) == (S(2)/7, [(x+S(1)/2, 1)])
+
+        if y is not None:
+
+            assert factor(x*y+1) == (S(1), [(x*y+1, 1)])
+            assert factor(x*y) == (S(1), [(x, 1), (y, 1)])
+
+            if characteristic_zero:
+                # primitive factors over Z for Z and Q.
+                assert factor(2*x + y) == (S(1), [(2*x + y, 1)])
+            else:
+                # monic factors over Z/pZ and GF(p^d)
+                assert factor(2*x + y) == (S(1)/2, [(x + y/2, 1)])
+
+
 def _all_matrices():
     """Return a list of matrix types and scalar types."""
     R163 = flint.fmpz_mod_ctx(163)
@@ -4152,6 +4231,8 @@ def test_all_tests():
     test_division_poly,
     test_division_matrix,
 
+    test_factor_poly_mpoly,
+
     test_polys,
     test_mpolys,
 

src/flint/types/fmpq_mpoly.pyx

@@ -940,7 +940,7 @@ cdef class fmpq_mpoly(flint_mpoly):
             c = fmpz.__new__(fmpz)
             fmpz_init_set((<fmpz>c).val, &fac.exp[i])
 
-            res[i] = (u, c)
+            res[i] = (u, int(c))
 
         c = fmpq.__new__(fmpq)
         fmpq_set((<fmpq>c).val, fac.constant)

src/flint/types/fmpq_poly.pyx

@@ -393,26 +393,41 @@ cdef class fmpq_poly(flint_poly):
         fmpq_poly_xgcd(res1.val, res2.val, res3.val, self.val, (<fmpq_poly>other).val)
         return (res1, res2, res3)
 
-    def factor(self):
+    def factor(self, *, monic=False):
         """
         Factors *self* into irreducible polynomials. Returns (*c*, *factors*)
         where *c* is the leading coefficient and *factors* is a list of
-        (*poly*, *exp*) pairs with all *poly* monic.
+        (*poly*, *exp*).
 
             >>> fmpq_poly.legendre_p(5).factor()
-            (63/8, [(x, 1), (x^4 + (-10/9)*x^2 + 5/21, 1)])
+            (1/8, [(x, 1), (63*x^4 + (-70)*x^2 + 15, 1)])
             >>> (fmpq_poly([1,-1],10) ** 5 * fmpq_poly([1,2,3],7)).factor()
+            (-1/700000, [(3*x^2 + 2*x + 1, 1), (x + (-1), 5)])
+
+        Since python-flint 0.7.0 this returns primitive denominator-free
+        factors consistent with ``fmpq_mpoly.factor()``. In previous versions
+        of python-flint all factors were made monic. Pass ``monic=True`` to get
+        monic factors instead.
+
+            >>> fmpq_poly.legendre_p(5).factor(monic=True)
+            (63/8, [(x, 1), (x^4 + (-10/9)*x^2 + 5/21, 1)])
+            >>> (fmpq_poly([1,-1],10) ** 5 * fmpq_poly([1,2,3],7)).factor(monic=True)
             (-3/700000, [(x^2 + 2/3*x + 1/3, 1), (x + (-1), 5)])
 
         """
         c, fac = self.numer().factor()
         c = fmpq(c)
-        for i in range(len(fac)):
-            base, exp = fac[i]
-            lead = base[base.degree()]
-            base = fmpq_poly(base, lead)
-            c *= lead ** exp
-            fac[i] = (base, exp)
+
+        if monic:
+            for i in range(len(fac)):
+                base, exp = fac[i]
+                lead = base[base.degree()]
+                base = fmpq_poly(base, lead)
+                c *= lead ** exp
+                fac[i] = (base, exp)
+        else:
+            fac = [(fmpq_poly(f), m) for f, m in fac]
+
         return c / self.denom(), fac
 
     def sqrt(self):

src/flint/types/fmpz_mpoly.pyx

@@ -919,7 +919,7 @@ cdef class fmpz_mpoly(flint_mpoly):
             c = fmpz.__new__(fmpz)
             fmpz_set((<fmpz>c).val, &fac.exp[i])
 
-            res[i] = (u, c)
+            res[i] = (u, int(c))
 
         c = fmpz.__new__(fmpz)
         fmpz_set((<fmpz>c).val, fac.constant)

src/flint/types/nmod_poly.pyx

@@ -235,10 +235,17 @@ cdef class nmod_poly(flint_poly):
             >>> f.leading_coefficient()
             133
         """
+        cdef nmod c
+
         d = self.degree()
         if d < 0:
             return 0
-        return nmod_poly_get_coeff_ui(self.val, d)
+
+        c = nmod.__new__(nmod)
+        c.mod = self.val.mod
+        c.val = nmod_poly_get_coeff_ui(self.val, d)
+
+        return c
 
     def inverse_series_trunc(self, slong n):
         """