Speed up roots when multiplicities=False for polynomials over finite fields

[c0rydoras]

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 roots of polynomials over finite fields is slower than it could be #35556

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

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes. Current tests still pass*
• I have updated the documentation and checked the documentation preview.

* one test fails because of ordering?

File "src/sage/rings/finite_rings/integer_mod_ring.py", line 1777, in sage.rings.finite_rings.integer_mod_ring.IntegerModRing_generic._roots_univariate_polynomial
Failed example:
    (x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots(multiplicities=False)
Expected:
    [19, 20, 21]
Got:
    [21, 20, 19]

⌛ Dependencies

[user202729]

In the worst case i.e. f is already a product of linear polynomials, how slow is the pow compared to the factorization?

[c0rydoras]

In the worst case i.e. f is already a product of linear polynomials, how slow is the pow compared to the factorization?

ehm, idrk, but my guess would be the difference would be more or less negligible in most cases

[user202729]

Do some benchmark then?

Also does this work when the modulus is more than 64 bits? I think normal pow (without modulus) sometimes fail in that case.

[c0rydoras]

Also does this work when the modulus is more than 64 bits?

yes, see the example in the description (also we don't do pow without a modulus)

[c0rydoras]

In the worst case i.e. f is already a product of linear polynomials, how slow is the pow compared to the factorization?

After some testing, the difference is < 1ms, should be fine?

[user202729]

okay. Do add a test that would have been (very) slow before the change and fast (let's say <1s) after the change. You may want to set_random_seed() to make the output & runtime deterministic.

Does this actually support polynomials over any finite fields (e.g. GF(5^2)), or just finite fields that is a quotient of ℤ?

[c0rydoras]

My current implementation only works for finite fields that are a quotient of ℤ

[TheBlupper]

I don't think that's true, exact same concept applies just raise to the actual order of the larger field.

sage: p = random_prime(2**64)
sage: F = GF((p, 5))
sage: R.<x> = F[]
sage: f = R.random_element(degree=500)
sage: %time f.roots()
CPU times: user 5.04 s, sys: 10.1 ms, total: 5.05 s
Wall time: 5.04 s
[(1696449494167646524*z5^4 + 8871908504779766200*z5^3 + 75123468697369816*z5^2 + 8425597716018871725*z5 + 3317261663788709534,
  1),
 (12557525520458311272*z5^4 + 1975064827054420205*z5^3 + 11479291003439526237*z5^2 + 11561078139694077531*z5 + 7844731532039912253,
  1)]
sage: %time (pow(x, F.order(), f)-x).gcd(f).roots()
CPU times: user 1.93 s, sys: 301 μs, total: 1.93 s
Wall time: 1.93 s
[(1696449494167646524*z5^4 + 8871908504779766200*z5^3 + 75123468697369816*z5^2 + 8425597716018871725*z5 + 3317261663788709534,
  1),
 (12557525520458311272*z5^4 + 1975064827054420205*z5^3 + 11479291003439526237*z5^2 + 11561078139694077531*z5 + 7844731532039912253,
  1)]

[c0rydoras]

My implementation only works for finite fields that are a quotient of ℤ, respectively thats the only time its actually called AFAICT (as it is right now)

[c0rydoras]

and also the .modulus() won't work for those

something like

R = f.parent()
x = R.gen()
p = R.base_ring().order()
g = pow(x, p, f) - x
g = f.gcd(g)
return g._roots_from_factorization(g.factor(), False)

could (i think?), work for all finite fields?

[TheBlupper]

You can use .characteristic() on finite fields to get the prime, but as I show in my snippet you really only need the order of the field. I definitively think this should be implemented for general finite fields, but as @user202729 pointed out some proper benchmarks for different sizes and shapes of polynomials would be good first. Could e.g be that this should be skipped for small enough polynomials.

[TheBlupper]

Did you investigate calling out to pari and using polrootsmod, as mentioned in #35556 ? From my limited testing it is quite fast, they probably already do this technique.

[c0rydoras]

Using the polynomial from the test I wrote:

sage: set_random_seed(31337)
sage: p = random_prime(2^128)
sage: R.<x> = GF(p)[]
sage: f = R.random_element(degree=5000)
sage: f.roots(multiplicities=False)

sage: %timeit pari.polrootsmod(f)
1.81 s ± 43.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
sage: %timeit (pow(x, F.order(), f)-x).gcd(f).roots()
388 ms ± 3.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
sage:

I also tried some other stuff, and while its faster than just .roots() (without this change), its not faster than (pow(x, F.order(), f)-x).gcd(f).roots(), not even close

[user202729]

this documentation is incorrectly formatted (actually this wouldn't show up in the documentation anyway, but you should format it correctly regardless)

Search for TESTS: in the code base for examples.

[user202729]

only add 4 spaces of indentation relative to the line above.

[c0rydoras]

i just adapted what i saw in src/sage/rings/finite_rings/integer_mod_ring.py, is there some sort of style guide for this?

[user202729]

it's just RST. See also sage coding_basics docs page.

[user202729]

only add 4 spaces of indentation relative to the line above.

this comment is still unaddressed.

[user202729]

wait for someone to approve workflows (after the change above) I suppose. But otherwise looks reasonable.

[user202729]

apart from https://github.com/sagemath/sage/pull/40494/files#r2248968501 , if tests pass (wait for someone to approve workflows…), should be good. remaining issues are negligible.

[c0rydoras]

Locally tests are happy except for:

File "src/sage/rings/finite_rings/integer_mod_ring.py", line 1777, in sage.rings.finite_rings.integer_mod_ring.IntegerModRing_generic._roots_univariate_polynomial
Failed example:
    (x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots(multiplicities=False)
Expected:
    [19, 20, 21]
Got:
    [21, 20, 19]

this one with an ordering issue

[user202729]

I suggest changing that test to

    sage: sorted((x^6 + x^5 + 9*x^4 + 20*x^3 + 3*x^2 + 18*x + 7).roots(multiplicities=False))
    [19, 20, 21]

[c0rydoras]

After rebasing, this no longer works, because now _roots_univariate_polynomial is called with algorithm=None, which results in a NotImplementedError, and then we continue in the roots function and end up with the polynomial being "factored"

edit: not sure when/where but I had algorithm is None instead of algorithm is not None 😅

[user202729]

I said algorithm is not None. Makes sense right?

[c0rydoras]

oopsie, fixed now 😅