Merge pull request #189 from oscarbenjamin/pr_fmpq_poly_factor_prim

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/flintlib/nmod_poly_factor.pxd

@@ -6,9 +6,9 @@ cdef extern from "flint/nmod_poly_factor.h":
 
     ctypedef struct nmod_poly_factor_struct:
         nmod_poly_struct *p
-        long *exp
-        long num
-        long alloc
+        slong *exp
+        slong num
+        slong alloc
     ctypedef nmod_poly_factor_struct nmod_poly_factor_t[1]
 
 # from here on is parsed

src/flint/test/test_all.py

@@ -486,3 +486,17 @@ cdef class fmpq(flint_scalar):
             return v
         else:
             raise OverflowError("fmpq_pow_fmpz(): exponent too large")
+
+    def sqrt(self):
+        """
+        Return exact rational square root of self or raise an error.
+
+            >>> fmpq(9, 4).sqrt()
+            3/2
+            >>> fmpq(8).sqrt()
+            Traceback (most recent call last):
+                ...
+            flint.utils.flint_exceptions.DomainError: not a square number
+
+        """
+        return fmpq(self.numer().sqrt(), self.denom().sqrt())

src/flint/types/fmpq.pyx

@@ -858,6 +858,25 @@ cdef class fmpq_mpoly(flint_mpoly):
         fmpq_mpoly_total_degree_fmpz((<fmpz> res).val, self.val, self.ctx.val)
         return res
 
+    def leading_coefficient(self):
+        """
+        Leading coefficient in the monomial ordering.
+
+            >>> from flint import Ordering
+            >>> ctx = fmpq_mpoly_ctx(2, Ordering.lex, ['x', 'y'])
+            >>> x, y = ctx.gens()
+            >>> p = 2*x*y + 3*x + 4*y**2 + 5
+            >>> p
+            2*x*y + 3*x + 4*y^2 + 5
+            >>> p.leading_coefficient()
+            2
+
+        """
+        if fmpq_mpoly_is_zero(self.val, self.ctx.val):
+            return fmpq(0)
+        else:
+            return self.coefficient(0)
+
     def repr(self):
         return f"{self.ctx}.from_dict({self.to_dict()})"
 
@@ -906,7 +925,7 @@ cdef class fmpq_mpoly(flint_mpoly):
         if fmpq_mpoly_sqrt(res.val, self.val, self.ctx.val):
             return res
         else:
-            raise ValueError("polynomial is not a perfect square")
+            raise DomainError("polynomial is not a perfect square")
 
     def factor(self):
         """
@@ -940,7 +959,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)
@@ -979,7 +998,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_mpoly.pyx

@@ -176,6 +176,24 @@ cdef class fmpq_poly(flint_poly):
     def is_one(self):
         return <bint>fmpq_poly_is_one(self.val)
 
+    def leading_coefficient(self):
+        """
+        Returns the leading coefficient of the polynomial.
+
+        >>> f = fmpq_poly([1, 2, 3])
+        >>> f
+        3*x^2 + 2*x + 1
+        >>> f.leading_coefficient()
+        3
+        """
+        cdef fmpq x
+        cdef slong d
+        d = fmpq_poly_degree(self.val)
+        x = fmpq.__new__(fmpq)
+        if d >= 0:
+            fmpq_poly_get_coeff_fmpq(x.val, self.val, d)
+        return x
+
     def __call__(self, other):
         t = any_as_fmpz(other)
         if t is not NotImplemented:
@@ -393,28 +411,64 @@ 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 factor_squarefree(self):
+        """
+        Factors *self* into square-free polynomials. Returns (*c*, *factors*)
+        where *c* is the leading coefficient and *factors* is a list of
+        (*poly*, *exp*).
+
+            >>> x = fmpq_poly([0, 1])
+            >>> p = x**2 * (x/2 - 1)**2 * (x + 1)**3
+            >>> p
+            1/4*x^7 + (-1/4)*x^6 + (-5/4)*x^5 + 1/4*x^4 + 2*x^3 + x^2
+            >>> p.factor_squarefree()
+            (1/4, [(x^2 + (-2)*x, 2), (x + 1, 3)])
+            >>> p.factor()
+            (1/4, [(x, 2), (x + (-2), 2), (x + 1, 3)])
+
+        """
+        c, fac = self.numer().factor_squarefree()
+        c = fmpq(c) / self.denom()
+        fac = [(fmpq_poly(f), m) for f, m in fac]
+        return c, fac
+
     def sqrt(self):
         """
         Return the exact square root of this polynomial or ``None``.

src/flint/types/fmpq_poly.pyx

@@ -883,29 +883,102 @@ cdef class fmpz(flint_scalar):
             return self.bit_length()
 
     def isqrt(self):
+        """
+        Return square root rounded down.
+
+            >>> fmpz(9).isqrt()
+            3
+            >>> fmpz(8).isqrt()
+            2
+
+        """
+        cdef fmpz v
+
+        if fmpz_sgn(self.val) < 0:
+            raise DomainError("integer square root of a negative number")
+
+        v = fmpz()
+        fmpz_sqrt(v.val, self.val)
+        return v
+
+    def sqrt(self):
+        """
+        Return exact integer square root of self or raise an error.
+
+            >>> fmpz(9).sqrt()
+            3
+            >>> fmpz(8).sqrt()
+            Traceback (most recent call last):
+                ...
+            flint.utils.flint_exceptions.DomainError: not a square number
+
+        """
         cdef fmpz v
+
         if fmpz_sgn(self.val) < 0:
-            raise ValueError("integer square root of a negative number")
+            raise DomainError("integer square root of a negative number")
+
         v = fmpz()
         fmpz_sqrt(v.val, self.val)
+
+        c = fmpz()
+        fmpz_mul(c.val, v.val, v.val)
+        if not fmpz_equal(c.val, self.val):
+            raise DomainError("not a square number")
+
         return v
 
     def sqrtrem(self):
+        """
+        Return the integer square root of self and remainder.
+
+            >>> fmpz(9).sqrtrem()
+            (3, 0)
+            >>> fmpz(8).sqrtrem()
+            (2, 4)
+            >>> c = fmpz(123456789012345678901234567890)
+            >>> u, v = c.sqrtrem()
+            >>> u ** 2 + v == c
+            True
+
+        """
         cdef fmpz u, v
+
         if fmpz_sgn(self.val) < 0:
-            raise ValueError("integer square root of a negative number")
+            raise DomainError("integer square root of a negative number")
+
         u = fmpz()
         v = fmpz()
         fmpz_sqrtrem(u.val, v.val, self.val)
+
         return u, v
 
     # warning: m should be prime!
-    def sqrtmod(self, m):
+    def sqrtmod(self, p):
+        """
+        Return modular square root of self modulo *p* or raise an error.
+
+            >>> fmpz(10).sqrtmod(13)
+            6
+            >>> (6**2) % 13
+            10
+            >>> fmpz(11).sqrtmod(13)
+            Traceback (most recent call last):
+                ...
+            flint.utils.flint_exceptions.DomainError: modular square root does not exist
+
+        The modulus *p* must be a prime number.
+        """
         cdef fmpz v
+
         v = fmpz()
-        m = fmpz(m)
-        if not fmpz_sqrtmod(v.val, self.val, (<fmpz>m).val):
-            raise ValueError("unable to compute modular square root")
+        if fmpz_is_zero(self.val):
+            return v
+
+        p = fmpz(p)
+        if not fmpz_sqrtmod(v.val, self.val, (<fmpz>p).val):
+            raise DomainError("modular square root does not exist")
+
         return v
 
     def root(self, long n):

src/flint/types/fmpz.pyx

@@ -10,6 +10,7 @@ from flint.flintlib.fmpz cimport(
     fmpz_is_probabprime,
     fmpz_mul,
     fmpz_invmod,
+    fmpz_sqrtmod,
     fmpz_divexact,
     fmpz_gcd,
     fmpz_is_one,
@@ -29,6 +30,9 @@ from flint.types.fmpz cimport(
 cimport cython
 cimport libc.stdlib
 
+from flint.utils.flint_exceptions import DomainError
+
+
 cdef class fmpz_mod_ctx:
     r"""
     Context object for creating :class:`~.fmpz_mod` initalised 
@@ -578,3 +582,34 @@ cdef class fmpz_mod(flint_scalar):
             )
 
         return res
+
+    def sqrt(self):
+        """
+        Return the square root of this ``fmpz_mod`` or raise an exception.
+
+            >>> ctx = fmpz_mod_ctx(13)
+            >>> s = ctx(10).sqrt()
+            >>> s
+            fmpz_mod(6, 13)
+            >>> s * s
+            fmpz_mod(10, 13)
+            >>> ctx(11).sqrt()
+            Traceback (most recent call last):
+                ...
+            flint.utils.flint_exceptions.DomainError: no square root exists for 11 mod 13
+
+        The modulus must be prime.
+
+        """
+        cdef fmpz_mod v
+
+        v = fmpz_mod.__new__(fmpz_mod)
+        v.ctx = self.ctx
+
+        if fmpz_is_zero(self.val):
+            return v
+
+        if not fmpz_sqrtmod(v.val, self.val, self.ctx.val.n):
+            raise DomainError("no square root exists for {} mod {}".format(self, self.ctx.modulus()))
+
+        return v