gh-38722: Roots of polynomials mod n
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->
Presently, Sage implements root-finding over Z/n in a way that requires
time exponential in e, where e is the largest exponent of a prime
divisor of n. For example, while Sage is capable of the following, it
takes 8 minutes and 40 seconds on the laptop I'm using:
```python
R.<x> = Zmod(2^30)[]
(x^2 - 1).roots(multiplicities=False)
```
This patch implements a better algorithm. Like the present one, it uses
the Chinese Remainder Theorem to reduce to the case of prime powers.
Unlike the present algorithm, it finds roots modulo prime powers using
Hensel lifting. When the root is non-singular mod p, it uses Newton
iteration and converges quadratically. When the root is singular,
lifting is done one power at a time, and the running time is bounded by
the product of e and the maximum number of roots modulo any p^d with d ≤
e. With this patch, the above code snippet takes 32 milliseconds; even
lifting modulo 2^1000 takes less than a second.
Depends on #38720.
URL: https://github.com/sagemath/sage/pull/38722
Reported by: Kyle Hofmann
Reviewer(s): Kyle Hofmann, Vincent Macri
