sagemathgh-37443: Fix bugs and regression in performance of `any_root()`
    
This pull request aims to fix some issues which came with the
refactoring of `any_root()` in
sagemath#37170.

I am afraid that apart from a little more cleaning up and verbose
comments, this "fix" really relies on `f.roots()` rather than
`f.any_root()` when roots over extension fields are computed.

The core problem currently is that the proper workflow should be the
following:

- Given a polynomial $f$, find an irreducible factor of smallest degree.
When `ring` is `None` and `degree` is `None`, only linear factors are
searched for. This code is unchanged and works fast.
- When `degree` is `None` and `ring` is not `None`, the old code used to
perform `f.change_ring(ring).any_root()` and look for a linear factor in
`ring`. The issue is when `ring` is a field with a large degree,
arithmetic in this ring is very slow and root finding is very slow.
- When `degree` is not `None` then an irreducible of degree `degree` is
searched for, $f$, and then a root is found in an extension ring. This
is either the user supplied `ring`, or it is in the extension
`self.base_ring().extension(degree)`. Again, as this extension could be
large, this can be slow.

Now, the reason that there's a regression in speed, as explained in
sagemath#37359, is that the method before
refactoring simply called `f.roots(ring)[0]` when working for these
extensions. As the root finding is always performed by a C library this
is much faster than the python code which we have for `any_root()` and
so the more "pure" method I wrote simply is slower.

The real issue holding this method back is that if we are in some field
$GF(p^k)$ and `ring` is some extension $GF(p^{nk})$ then we are not
guaranteed a coercion between these rings with the current coercion
system. Furthermore, if we extend the base to some $GF(p^{dk})$ with
$dk$ dividing $nk$ we also have no certainty of a coercion from this
minimal extension to the user supplied ring.

As a result, a "good" fix is delayed until we figure out finite field
coercion. **Issue to come**.

Because of all of this, this PR goes back to the old behaviour, while
keeping the refactored and easier to read code. I hope one day in the
future to fix this.

Fixes sagemath#37359
Fixes sagemath#37442
Fixes sagemath#37445
Fixes sagemath#37471

**EDIT**

I have labeled this as critical as the slow down in the most extreme
cases causes code to complete hang when looking for a root. See for
example: sagemath#37442. I do not want
sagemath#37170 to be merged into 10.3
without this patch.

Additionally, a typo meant that a bug was introduced in the CZ
splitting. This has also been fixed.

Also, also in the refactoring the function behaved as intended (rather
than the strange behaviour from before) which exposed an issue
sagemath#37471 which I have fixed.
    
URL: sagemath#37443
Reported by: Giacomo Pope
Reviewer(s): Giacomo Pope, grhkm21, Lorenz Panny, Travis Scrimshaw

src/sage/rings/finite_rings/conway_polynomials.py

@@ -239,7 +239,7 @@ def polynomial(self, n):
         # TODO: something like the following
         # gcds = [n.gcd(d) for d in self.nodes.keys()]
         # xi = { m: (...) for m in gcds }
-        xi = {q: self.polynomial(n//q).any_root(K, n//q, assume_squarefree=True, assume_distinct_deg=True)
+        xi = {q: self.polynomial(n//q).any_root(K, n//q, assume_squarefree=True, assume_equal_deg=True)
               for q in n.prime_divisors()}
 
         # The following is needed to ensure that in the concrete instantiation

src/sage/rings/polynomial/polynomial_element.pyx

@@ -134,7 +134,6 @@ from sage.categories.morphism cimport Morphism
 from sage.misc.superseded import deprecation_cython as deprecation, deprecated_function_alias
 from sage.misc.cachefunc import cached_method
 
-
 cpdef is_Polynomial(f) noexcept:
     """
     Return ``True`` if ``f`` is of type univariate polynomial.
@@ -2151,6 +2150,18 @@ cdef class Polynomial(CommutativePolynomial):
             True
             sage: f % factor == 0
             True
+
+        TESTS:
+
+        Ensure that :issue:`37445` is fixed::
+
+            sage: R.<x> = GF(13)[]
+            sage: def irr(d, R): return f.monic() if (f := R.random_element(d)).is_irreducible() else irr(d, R)
+            sage: f = prod(irr(6, R) for _ in range(10))
+            sage: irr = f._cantor_zassenhaus_split_to_irreducible(6)
+            sage: assert irr.degree() == 6
+            sage: assert f % irr == 0
+            sage: assert irr.is_irreducible()
         """
         R = self.parent()
         q = self.base_ring().order()
@@ -2160,8 +2171,8 @@ cdef class Polynomial(CommutativePolynomial):
             return self
 
         # We expect to succeed with greater than 1/2 probability,
-        # so if we try 1000 times and fail, there's a bug somewhere.
-        for _ in range(1000):
+        # so if we try 100 times and fail, there's a bug somewhere.
+        for _ in range(100):
             # Sample a polynomial "uniformly" from R
             # TODO: once #37118 has been merged, this can be made cleaner,
             #       as we will actually have access to uniform sampling.
@@ -2174,7 +2185,7 @@ cdef class Polynomial(CommutativePolynomial):
 
             # Need to handle odd and even characteristic separately
             if q % 2:
-                h = self.gcd(pow(T, (q-1)//2, self) - 1)
+                h = self.gcd(pow(T, (q**degree-1)//2, self) - 1)
             else:
                 # Compute the trace of T with field of order 2^k
                 # sum T^(2^i) for i in range (degree * k)
@@ -2200,20 +2211,23 @@ cdef class Polynomial(CommutativePolynomial):
         # If you are reaching this error, chances are there's a bug in the code.
         raise AssertionError(f"no splitting of degree {degree} found for {self}")
 
-    def _any_irreducible_factor_squarefree(self, degree=None):
+    def _any_irreducible_factor_squarefree(self, degree=None, ext_degree=None):
         """
-        Helper function for any_irreducible_factor which computes
+        Helper function for :meth:`any_irreducible_factor` which computes
         an irreducible factor from self, assuming the input is
         squarefree.
 
         Does this by first computing the distinct degree factorisations
         of self and then finds a factor with Cantor-Zassenhaus
         splitting.
 
-        If degree is not None, then only irreducible factors of degree
-        `degree` are searched for, otherwise the smallest degree factor
+        If ``degree`` is not ``None``, then only irreducible factors of degree
+        ``degree`` are searched for, otherwise the smallest degree factor
         is found.
 
+        If ``ext_degree`` is not ``None``, then only irreducible factors whose
+        degree divides ``ext_degree`` are returned.
+
         EXAMPLES::
 
             sage: # needs sage.rings.finite_rings
@@ -2262,10 +2276,15 @@ cdef class Polynomial(CommutativePolynomial):
 
         # Otherwise we use the smallest possible d value
         for (poly, d) in self._distinct_degree_factorisation_squarefree():
-            return poly._cantor_zassenhaus_split_to_irreducible(d)
-        raise ValueError(f"no irreducible factor could be computed from {self}")
+            if ext_degree is None:
+                return poly._cantor_zassenhaus_split_to_irreducible(d)
+            elif ZZ(d).divides(ext_degree):
+                return poly._cantor_zassenhaus_split_to_irreducible(d)
+            if d > ext_degree:
+                raise ValueError(f"no irreducible factor of degree {degree} dividing {ext_degree} could be computed from {self}")
+        raise AssertionError(f"no irreducible factor could be computed from {self}")
 
-    def any_irreducible_factor(self, degree=None, assume_squarefree=False, assume_distinct_deg=False):
+    def any_irreducible_factor(self, degree=None, assume_squarefree=False, assume_equal_deg=False, ext_degree=None):
         """
         Return an irreducible factor of this polynomial.
 
@@ -2281,11 +2300,15 @@ cdef class Polynomial(CommutativePolynomial):
           Used for polynomials over finite fields.  If ``True``,
           this polynomial is assumed to be squarefree.
 
-        - ``assume_distinct_deg`` (boolean) -- (default: ``False``).
+        - ``assume_equal_deg`` (boolean) -- (default: ``False``).
           Used for polynomials over finite fields.  If ``True``,
           this polynomial is assumed to be the product of irreducible
           polynomials of degree equal to ``degree``.
 
+        - ``ext_degree`` -- positive integer or ``None`` (default);
+          used for polynomials over finite fields. If not ``None`` only returns
+          irreducible factors of ``self`` whose degree divides ``ext_degree``.
+
         EXAMPLES::
 
             sage: # needs sage.rings.finite_rings
@@ -2328,9 +2351,9 @@ cdef class Polynomial(CommutativePolynomial):
             sage: F = GF(163)
             sage: R.<x> = F[]
             sage: h = (x + 57) * (x + 98) * (x + 117) * (x + 145)
-            sage: h.any_irreducible_factor(degree=1, assume_distinct_deg=True)   # random
+            sage: h.any_irreducible_factor(degree=1, assume_equal_deg=True)   # random
             x + 98
-            sage: h.any_irreducible_factor(assume_distinct_deg=True)
+            sage: h.any_irreducible_factor(assume_equal_deg=True)
             Traceback (most recent call last):
             ...
             ValueError: degree must be known if distinct degree factorisation is assumed
@@ -2359,7 +2382,7 @@ cdef class Polynomial(CommutativePolynomial):
             if degree < 1:
                 raise ValueError(f"{degree = } must be positive")
 
-        if assume_distinct_deg and degree is None:
+        if assume_equal_deg and degree is None:
             raise ValueError("degree must be known if distinct degree factorisation is assumed")
 
         # When not working over a finite field, do the simple thing of factoring.
@@ -2397,27 +2420,30 @@ cdef class Polynomial(CommutativePolynomial):
 
         # If we know the polynomial is square-free, we can start here
         if assume_squarefree:
-            if assume_distinct_deg:
+            if assume_equal_deg:
                 return self._cantor_zassenhaus_split_to_irreducible(degree)
-            return self._any_irreducible_factor_squarefree(degree)
+            return self._any_irreducible_factor_squarefree(degree, ext_degree)
 
         # Otherwise we compute the squarefree decomposition and check each
         # polynomial for a root. If no poly has a root, we raise an error.
         SFD = self.squarefree_decomposition()
         SFD.sort()
         for poly, _ in SFD:
             try:
-                return poly._any_irreducible_factor_squarefree(degree)
+                return poly._any_irreducible_factor_squarefree(degree, ext_degree)
             except ValueError:
                 pass
 
         # If degree has been set, there could just be no factor of the desired degree
         if degree:
             raise ValueError(f"polynomial {self} has no irreducible factor of degree {degree}")
+        # If ext_degree has been set, then there may be no irreducible factor of degree dividing ext_degree
+        if ext_degree:
+            raise ValueError(f"polynomial {self} has no irreducible factor of degree dividing {ext_degree}")
         # But if any degree is allowed then there should certainly be a factor if self has degree > 0
         raise AssertionError(f"no irreducible factor was computed for {self}. Bug.")
 
-    def any_root(self, ring=None, degree=None, assume_squarefree=False, assume_distinct_deg=False):
+    def any_root(self, ring=None, degree=None, assume_squarefree=False, assume_equal_deg=False):
         """
         Return a root of this polynomial in the given ring.
 
@@ -2437,14 +2463,26 @@ cdef class Polynomial(CommutativePolynomial):
           finite fields.  If ``True``, this polynomial is assumed to be
           squarefree.
 
-        - ``assume_distinct_deg`` (bool) -- Used for polynomials over
+        - ``assume_equal_deg`` (bool) -- Used for polynomials over
           finite fields.  If ``True``, all factors of this polynomial
-          are assumed to have degree ``degree``.
+          are assumed to have degree ``degree``. Note that ``degree``
+          must be set.
 
         .. WARNING::
 
             Negative degree input will be deprecated. Instead use
-            ``assume_distinct_deg``.
+            ``assume_equal_deg``.
+
+        .. NOTE::
+
+            For finite fields, ``any_root()`` is non-deterministic when
+            finding linear roots of a polynomial over the base ring.
+            However, if ``degree`` is greater than one, or ``ring`` is an
+            extension of the base ring, then the root computed is found
+            by attempting to return a root after factorisation. Roots found
+            in this way are deterministic. This may change in the future.
+            For all other rings or fields, roots are found by first
+            fully-factoring ``self`` and the output is deterministic.
 
         EXAMPLES::
 
@@ -2554,38 +2592,76 @@ cdef class Polynomial(CommutativePolynomial):
         # When not working over a finite field, do the simple thing of factoring for
         # roots and picking the first root. If none are available, raise an error.
         from sage.categories.finite_fields import FiniteFields
-        if not self.base_ring() in FiniteFields():
-            rs = self.roots(ring=ring, multiplicities=False)
-            if rs:
-                return rs[0]
-            raise ValueError(f"polynomial {self} has no roots")
+        if self.base_ring() not in FiniteFields():
+            if ring not in FiniteFields():
+                rs = self.roots(ring=ring, multiplicities=False)
+                if rs:
+                    return rs[0]
+                raise ValueError(f"polynomial {self} has no roots")
+
+        # Ensure that a provided ring is appropriate for the function. From the
+        # above we know it is either None or a finite field. When it's a finite
+        # field we ensure there's a coercion from the base ring to ring.
+        if ring is not None:
+            if ring.coerce_map_from(self.base_ring()) is None:
+                raise ValueError(f"no coercion map can be computed from {self.base_ring()} to {ring}")
 
         # When the degree is none, we only look for a linear factor
         if degree is None:
-            # if a ring is given try and coerce the polynomial into this ring
-            if ring is not None:
+            # When ring is None, we attempt to find a linear factor of self
+            if ring is None:
                 try:
-                    self = self.change_ring(ring)
+                    f = self.any_irreducible_factor(degree=1, assume_squarefree=assume_squarefree)
                 except ValueError:
-                    raise(f"cannot coerce polynomial {self} to the new ring: {ring}")
+                    raise ValueError(f"no root of polynomial {self} can be computed")
+                return - f[0] / f[1]
 
-            # try and find a linear irreducible polynomial from f to compute a root
+            # When we have a ring, then we can find an irreducible factor of degree `d` providing
+            # that d divides the degree of the extension from the base ring to the given ring
+            allowed_extension_degree = ring.degree() // self.base_ring().degree()
             try:
-                f = self.any_irreducible_factor(degree=1, assume_squarefree=assume_squarefree)
+                f = self.any_irreducible_factor(assume_squarefree=assume_squarefree, ext_degree=allowed_extension_degree)
             except ValueError:
-                raise ValueError(f"no root of polynomial {self} can be computed")
-
-            return - f[0] / f[1]
+                raise ValueError(f"no root of polynomial {self} can be computed over the ring {ring}")
+            # When d != 1 we then find the smallest extension
+            # TODO: What we should do here is compute some minimal
+            #       extension F_ext = self.base_ring().extension(d, names="a") and find a
+            #       root here and then coerce this root into the parent ring. This means we
+            #       would work with the smallest possible extension.
+            #       However, if we have some element of GF(p^k) and we try and coerce this to
+            #       some element GF(p^(k*n)) this can fail, even though mathematically it
+            #       should be fine.
+            # TODO: Additionally, if the above was solved, it would be faster to extend the base
+            #       ring with the irreducible factor however, if the base ring is an extension
+            #       then the type of self.base_ring().extension(f) is a Univariate Quotient Polynomial Ring
+            #       and not a finite field.
+
+            # When f has degree one we simply return the roots
+            # TODO: should we write something fast for degree two using
+            #       the quadratic formula?
+            if f.degree().is_one():
+                root = - f[0] / f[1]
+                return ring(root)
+
+            # TODO: The proper thing to do here would be to call
+            #       return f.change_ring(ring).any_root()
+            #       but as we cannot work in the minimal extension (see above) working
+            #       in the extension for f.change_ring(ring).any_root() is almost always
+            #       much much much slower than using the call for roots() which uses
+            #       C library bindings for all finite fields.
+            #       Until the coercion system for finite fields works better,
+            #       this will be the most performant
+            return f.roots(ring, multiplicities=False)[0]
 
         # The old version of `any_root()` allowed degree < 0 to indicate that the input polynomial
         # had a distinct degree factorisation, we pass this to any_irreducible_factor as a bool and
         # ensure that the degree is positive.
         degree = ZZ(degree)
         if degree < 0:
             from sage.misc.superseded import deprecation
-            deprecation(37170, "negative ``degree`` will be disallowed. Instead use the bool `assume_distinct_deg`.")
+            deprecation(37170, "negative ``degree`` will be disallowed. Instead use the bool `assume_equal_deg`.")
             degree = -degree
-            assume_distinct_deg = True
+            assume_equal_deg = True
 
         # If a certain degree is requested, then we find an irreducible factor of degree `degree`
         # use this to compute a field extension and return the generator as root of this polynomial
@@ -2594,7 +2670,7 @@ cdef class Polynomial(CommutativePolynomial):
         try:
             f = self.any_irreducible_factor(degree=degree,
                                             assume_squarefree=assume_squarefree,
-                                            assume_distinct_deg=assume_distinct_deg)
+                                            assume_equal_deg=assume_equal_deg)
         except ValueError:
             raise ValueError(f"no irreducible factor of degree {degree} can be computed from {self}")
 
@@ -2617,13 +2693,17 @@ cdef class Polynomial(CommutativePolynomial):
             #       FiniteField type if the base field is a non-prime field,
             #       so this slower option is chosen to ensure the root is
             #       over explicitly a FiniteField type.
-            ring = self.base_ring().extension(f.degree(), names="a")
-
-        # Now we look for a linear root of this irreducible polynomial of degree `degree`
-        # over the user supplied ring or the extension we just computed. If the user sent
-        # a bad ring here of course there may be no root found.
-        f = f.change_ring(ring)
-        return f.any_root()
+            ring = self.base_ring().extension(degree, names="a")
+
+        # TODO: The proper thing to do here would be to call
+        #       return f.change_ring(ring).any_root()
+        #       but as we cannot work in the minimal extension (see above) working
+        #       in the extension for f.change_ring(ring).any_root() is almost always
+        #       much much much slower than using the call for roots() which uses
+        #       C library bindings for all finite fields.
+        #       Until the coercion system for finite fields works better,
+        #       this will be the most performant
+        return f.roots(ring, multiplicities=False)[0]
 
     def __truediv__(left, right):
         r"""

src/sage/rings/polynomial/polynomial_quotient_ring.py

@@ -2098,7 +2098,7 @@ def _isomorphic_ring(self):
             base_image = self.base_ring().modulus().change_ring(isomorphic_ring).any_root()
             base_to_isomorphic_ring = self.base_ring().hom([isomorphic_ring(base_image)])
             modulus = self.modulus().map_coefficients(base_to_isomorphic_ring)
-            gen = modulus.any_root(assume_squarefree=True, degree=1, assume_distinct_deg=True)
+            gen = modulus.any_root(assume_squarefree=True, degree=1, assume_equal_deg=True)
 
             homspace = Hom(self, isomorphic_ring)
             to_isomorphic_ring = homspace.__make_element_class__(SetMorphism)(homspace,

src/sage/schemes/elliptic_curves/cm.py

@@ -233,6 +233,13 @@ def is_HCP(f, check_monic_irreducible=True):
         sage: all(not is_HCP(hilbert_class_polynomial(D) + 1)
         ....:     for D in srange(-4,-100,-1) if D.is_discriminant())
         True
+
+    Ensure that :issue:`37471` is fixed::
+
+        sage: from sage.schemes.elliptic_curves.cm import is_HCP
+        sage: set_random_seed(297388353221545796156853787333338705098)
+        sage: is_HCP(hilbert_class_polynomial(-55))
+        -55
     """
     zero = ZZ(0)
     # optional check that input is monic and irreducible
@@ -264,14 +271,15 @@ def is_HCP(f, check_monic_irreducible=True):
         # Compute X^p-X mod fp
         z = fp.parent().gen()
         r = pow(z, p, fp) - z
-        d = r.gcd(fp).degree()  # number of roots mod p
+        r = r.gcd(fp)
+        d = r.degree()  # number of roots mod p
         if d == 0:
             continue
-        if not fp.is_squarefree():
+        if not r.is_squarefree():
             continue
         if d < h and d not in h2list:
             return zero
-        jp = fp.any_root(degree=1, assume_squarefree=True, assume_distinct_deg=True)
+        jp = r.any_root(degree=1, assume_squarefree=True, assume_equal_deg=True)
         E = EllipticCurve(j=jp)
         if E.is_supersingular():
             continue