gh-40494: Speed up `roots` when `multiplicities=False` for polynomials over finite fields
    
Speed up `f.roots` when `multiplicities=False` for univariate
polynomials over finite fields

- as seen in this blog post: https://adib.au/2025/lance-hard/#speedup-
by-using-gcd
- also mentioned here #35556

```sage
sage: p = random_prime(2^200)
sage: F = GF(p)
sage: R.<x> = F[]
sage: f = R.random_element(degree=20000)
sage: %time f.roots(multiplicities=False) # would hang / be very slow
before this change
CPU times: user 5.24 s, sys: 16.4 ms, total: 5.26 s
Wall time: 5.26 s
[205828247017582431219769663595709973936079117296023914605767] # just
some roots, will ofc be different every run :)
```


### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation and checked the documentation
preview.
    
URL: #40494
Reported by: Arthur
Reviewer(s): Arthur, user202729

Author: Release Manager

src/sage/rings/finite_rings/finite_field_base.pyx

@@ -2120,6 +2120,83 @@ cdef class FiniteField(Field):
         python_int = int.from_bytes(input_bytes, byteorder=byteorder)
         return self.from_integer(python_int)
 
+    def _roots_univariate_polynomial(self, f, ring, multiplicities, algorithm=None):
+        r"""
+        Return the roots of the univariate polynomial ``f``.
+
+        INPUT:
+
+        - ``f`` -- a polynomial defined over this field
+
+        - ``ring`` -- the ring to find roots in.
+
+        - ``multiplicities`` -- bool (default: ``True``). If ``True``, return
+          list of pairs `(r, n)`, where `r` is a root and `n` is its
+          multiplicity. If ``False``, just return the unique roots, with no
+          information about multiplicities.
+
+        - ``algorithm`` -- ignored
+
+        TESTS:
+
+        We can take the roots of a polynomial defined over a finite field::
+
+            sage: set_random_seed(31337)
+            sage: p = random_prime(2^128)
+            sage: R.<x> = Zmod(p)[]
+            sage: f = R.random_element(degree=15)
+            sage: f.roots()
+            [(117558869610275297997958296126212805270, 1)]
+
+        We can take the roots of a polynomial defined over a finite field without multiplicities::
+
+            sage: set_random_seed(31337)
+            sage: p = random_prime(2^128)
+            sage: R.<x> = Zmod(p)[]
+            sage: f = R.random_element(degree=150)
+            sage: f.roots(multiplicities=False)
+            [116560079209701720510648792531840294827]
+
+        We can take the roots of a polynomial defined over a finite field extension::
+
+            sage: set_random_seed(31337)
+            sage: F.<a> = GF((2, 10))
+            sage: R.<x> = F[]
+            sage: f = R.random_element(degree=10)
+            sage: f.roots()
+            [(a^9 + a^8 + a^6 + a^4 + a^2, 1)]
+            sage: f.roots(multiplicities=False)
+            [a^9 + a^8 + a^6 + a^4 + a^2]
+
+        We can take the roots of a high degree polynomial in a reasonable time::
+
+            sage: set_random_seed(31337)
+            sage: p = random_prime(2^128)
+            sage: F = GF(p)
+            sage: R.<x> = F[]
+            sage: f = R.random_element(degree=10000)
+            sage: f.roots(multiplicities=False)
+            [65940671326230628578511607550463701471]
+        """
+        if algorithm is not None:
+            raise NotImplementedError
+
+        K = f.base_ring()
+        L = K if ring is None else ring
+
+        if K != L:
+            raise NotImplementedError
+
+        if multiplicities:
+            return f._roots_from_factorization(f.factor(), multiplicities)
+        else:
+            R = f.parent()
+            x = R.gen()
+            p = K.order()
+            g = pow(x, p, f) - x
+            g = f.gcd(g)
+            return g._roots_from_factorization(g.factor(), False)
+
 
 def unpickle_FiniteField_ext(_type, order, variable_name, modulus, kwargs):
     r"""

src/sage/rings/finite_rings/integer_mod_ring.py

@@ -1774,7 +1774,7 @@ def _roots_univariate_polynomial(self, f, ring=None, multiplicities=True, algori
             [0, 32, 9]
             sage: (x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots()
             [(19, 1), (20, 2), (21, 3)]
-            sage: (x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots(multiplicities=False)
+            sage: sorted((x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots(multiplicities=False))
             [19, 20, 21]
 
         We can find roots without multiplicities over a ring whose modulus is
@@ -1911,6 +1911,15 @@ def _roots_univariate_polynomial(self, f, ring=None, multiplicities=True, algori
             sage: R.<x> = Zmod(100)[]
             sage: (x^2 - 1).roots(Zmod(99), multiplicities=False) == (x^2 - 1).change_ring(Zmod(99)).roots(multiplicities=False)
             True
+
+        We can find roots of high degree polynomials in a reasonable time:
+
+            sage: set_random_seed(31337)
+            sage: p = random_prime(2^128)
+            sage: R.<x> = Zmod(p)[]
+            sage: f = R.random_element(degree=5000)
+            sage: f.roots(multiplicities=False)
+            [107295314027801680550847462044796892009, 75545907600948005385964943744536832524]
         """
 
         # This function only supports roots in an IntegerModRing
@@ -1926,7 +1935,7 @@ def _roots_univariate_polynomial(self, f, ring=None, multiplicities=True, algori
                     " implemented (try the multiplicities=False option)"
                 )
             # Roots of non-zero polynomial over finite fields by factorization
-            return f._roots_from_factorization(f.factor(), multiplicities)
+            return f.change_ring(f.base_ring().field()).roots(multiplicities=multiplicities)
 
         # Zero polynomial is a base case
         if deg < 0:
@@ -1935,7 +1944,12 @@ def _roots_univariate_polynomial(self, f, ring=None, multiplicities=True, algori
 
         # Finite fields are a base case
         if self.is_field():
-            return f._roots_from_factorization(f.factor(), False)
+            return list(
+                map(
+                    f.base_ring(),
+                    f.change_ring(f.base_ring().field()).roots(multiplicities=False),
+                )
+            )
 
         # Otherwise, find roots modulo each prime power
         fac = self.factored_order()

src/sage/rings/polynomial/polynomial_element.pyx

@@ -8377,8 +8377,8 @@ cdef class Polynomial(CommutativePolynomial):
             sage: A = PolynomialRing(R, 'y')
             sage: y = A.gen()
             sage: f = 10*y^2 - y^3 - 9
-            sage: f.roots(multiplicities=False)                                         # needs sage.libs.pari
-            [1, 0, 3, 6]
+            sage: sorted(f.roots(multiplicities=False))                                 # needs sage.libs.pari
+            [0, 1, 3, 6]
 
         An example over the complex double field (where root finding is
         fast, thanks to NumPy)::