Slow root finding over GF(2^k)[] quotient rings

[grhkm21]

Steps To Reproduce

sage: K.<z16> = GF(2^16)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: 
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: poly = R_ext.random_monic_element(2)
....: print("poly:", poly)
....: 
....: poly.roots()

Expected Behavior

It should take a short time, especially since here u is irreducible, so K_ext is in theory isomorphic to GF(2^32), and over that field it takes no time to find roots, just do hensel lift or whatever:

sage: K.<x> = GF(2^32)[]
....: %timeit K.random_monic_element(2).roots()
1.37 ms ± 55.2 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Actual Behavior

It's slow... Even changing the base field from 2^16 to 2^8 takes forever. Moreover, the time is spent on this weird self._roots_from_factorization(self.factor(), multiplicities) call. I suspect it is because I specify the modulus in the K_ext construction, so it has to compute some field isomorphism. Still, I don't see how 8-bit computations will take this long.

For more details, here's the timings:

sage: K.<z16> = GF(2^8)
....: R.<x> = K[]
....: u = R.irreducible_element(2)
....: 
....: K_ext = K.extension(modulus=u, names="a")
....: R_ext.<y_ext> = K_ext[]
....: %time R_ext.random_monic_element(2).roots()
....: %time _ = [R_ext.random_monic_element(2).roots() for _ in range(10^3)]
CPU times: user 43.8 s, sys: 163 ms, total: 44 s
Wall time: 44.4 s
[]
CPU times: user 3.91 s, sys: 15 ms, total: 3.93 s
Wall time: 3.94 s

As seen, a ridiculous amount of time is spent on some kind of field isomorphism / data initialisation.

Additional Information

I tried tracing the calls myself, but as mentioned, it goes into some weird call path _roots_from_factorization that I can't understand. Any help with either debugging the above or redirect it to use GF(2^32) like I suggested, would be appreciated.

For more additional information, this behaviour only seem to occur for $K = \mathbb{F}_q$ where $q = 2^k \geq 2^8$. Also once it finishes (which takes several tens of seconds), the root finding itself is very fast. To observe this, run while True: R_ext.random_monic_element(2).roots(). This behaviour also doesn't happen, even for weird fields like $K = \mathbb{F}_{103^{13}}$ the whole process takes < 1s.

[GiacomoPope]

The core of this issue is that when the isomorphism is computed, the function _isomorphic_ring() (which is cached) requires finding some roots of the modulus of the base ring:

            # the map to GF(N) maps our generator to a root of our modulus in the isomorphic_ring
            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)

Running %prun -s funtime f.roots() as you have above, we see that any_root() is what is taking up all the time:

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    5.152    5.152 {built-in method builtins.exec}
        1    0.000    0.000    5.152    5.152 <string>:1(<module>)
        1    0.000    0.000    5.152    5.152 {method 'roots' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
        1    0.000    0.000    5.152    5.152 polynomial_quotient_ring.py:1936(_factor_univariate_polynomial)
        1    0.003    0.003    5.143    5.143 polynomial_quotient_ring.py:1979(_isomorphic_ring)
        2    1.578    0.789    5.081    2.540 {method 'any_root' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
     7514    0.556    0.000    3.300    0.000 polynomial_ring.py:1363(random_element)

Looking at the any_root() function itself shows some amusing timings:

sage: F = GF(2^16)
sage: R.<x> = F[]
sage: f = R.random_element(degree=4)
sage:
sage: %time f.roots()
CPU times: user 1.91 ms, sys: 10 µs, total: 1.92 ms
Wall time: 1.94 ms
[(z16^14 + z16^12 + z16^10 + z16^9 + z16^7 + z16^6 + z16^5 + z16^4 + z16^2 + z16,
  1),
 (z16^15 + z16^14 + z16^13 + z16^12 + z16^11 + z16^9 + z16^8 + z16^2, 1)]
sage:
sage: %time f.any_root()
CPU times: user 6.52 s, sys: 11.9 ms, total: 6.53 s
Wall time: 6.54 s
z16^14 + z16^12 + z16^10 + z16^9 + z16^7 + z16^6 + z16^5 + z16^4 + z16^2 + z16

So we can look at what happens for odd characteristic:

sage: F = GF(3^16)
sage: R.<x> = F[]
sage: f = R.random_element(degree=4)
sage: %time f.roots()
CPU times: user 6.71 ms, sys: 1.56 ms, total: 8.27 ms
Wall time: 12.5 ms
[(2*z16^15 + z16^14 + 2*z16^13 + 2*z16^12 + 2*z16^11 + z16^10 + 2*z16^9 + z16^8 + 2*z16^7 + 2*z16^6 + z16^5 + 2*z16^4 + z16^2 + 2*z16,
  1)]
sage: %time f.any_root()
CPU times: user 9.82 ms, sys: 357 µs, total: 10.2 ms
Wall time: 9.87 ms
2*z16^15 + z16^14 + 2*z16^13 + 2*z16^12 + 2*z16^11 + z16^10 + 2*z16^9 + z16^8 + 2*z16^7 + 2*z16^6 + z16^5 + 2*z16^4 + z16^2 + 2*z16

and we see that this issue really seems to be on even characteristic only...

Here is the current implementation of any_root(), snippet to show the even / odd char parts:

    def any_root(self, ring=None, degree=None, assume_squarefree=False):
       #      
       # SNIPPED
       #
            q = self.base_ring().order()
            if q % 2 == 0:
                while True:
                    T = R.random_element(2*degree-1)
                    if T == 0:
                        continue
                    T = T.monic()
                    C = T
                    for i in range(degree-1):
                        C = T + pow(C,q,self)
                    h = self.gcd(C)
                    hd = h.degree()
                    if hd != 0 and hd != self.degree():
                        if 2*hd <= self.degree():
                            return h.any_root(ring, -degree, True)
                        else:
                            return (self//h).any_root(ring, -degree, True)
            else:
                while True:
                    T = R.random_element(2*degree-1)
                    if T == 0:
                        continue
                    T = T.monic()
                    h = self.gcd(pow(T, Integer((q**degree-1)/2), self)-1)
                    hd = h.degree()
                    if hd != 0 and hd != self.degree():
                        if 2*hd <= self.degree():
                            return h.any_root(ring, -degree, True)
                        else:
                            return (self//h).any_root(ring, -degree, True)

This raises the question, should any_root be fixed, or should we simply use choice(factor(f)) with suitable weighting, as most of the time (all the time??) we have C library support for the factoring.

This slowness seems to be related to the issues #21998 and #16162