Trac #34512: compute list CRT via tree

Currently, `CRT_list()` works by ''folding'' the input from one side. In
typical cases of interest, it is much more efficient to use a binary-
tree structure instead. (This is similar to how `prod()` is
implemented.)

Example:
{{{#!sage
sage: ms = list(primes(10^5))
sage: xs = [randrange(m) for m in ms]
sage: %timeit CRT_list(xs, ms)
}}}

Sage 9.7.rc0:
{{{
7.42 s ± 20 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
}}}

This branch:
{{{
86.5 ms ± 169 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
}}}

Similar speedups can be observed for polynomials; example (a length‑1024
inverse DFT):
{{{
sage: F = GF(65537)
sage: a = F(1111)
sage: assert a.multiplicative_order() == 1024
sage: R.<x> = F[]
sage: ms = [x - a^i for i in range(1024)]
sage: zs = [R(F.random_element()) for _ in ms]
sage: %timeit CRT_list(zs, ms)
}}}

Sage 9.7.rc0:
{{{
347 ms ± 863 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
}}}

This branch:
{{{
29.3 ms ± 37.1 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
}}}

URL: https://trac.sagemath.org/34512
Reported by: lorenz
Ticket author(s): Lorenz Panny
Reviewer(s): Kwankyu Lee

src/sage/arith/misc.py

@@ -3224,10 +3224,10 @@ def crt(a, b, m=None, n=None):
 CRT = crt
 
 
-def CRT_list(v, moduli):
-    r""" Given a list ``v`` of elements and a list of corresponding
+def CRT_list(values, moduli):
+    r""" Given a list ``values`` of elements and a list of corresponding
     ``moduli``, find a single element that reduces to each element of
-    ``v`` modulo the corresponding moduli.
+    ``values`` modulo the corresponding moduli.
 
     .. SEEALSO::
 
@@ -3289,22 +3289,37 @@ def CRT_list(v, moduli):
         sage: from gmpy2 import mpz
         sage: CRT_list([mpz(2),mpz(3),mpz(2)], [mpz(3),mpz(5),mpz(7)])
         23
+
+    Make sure we are not mutating the input lists::
+
+        sage: xs = [1,2,3]
+        sage: ms = [5,7,9]
+        sage: CRT_list(xs, ms)
+        156
+        sage: xs
+        [1, 2, 3]
+        sage: ms
+        [5, 7, 9]
     """
-    if not isinstance(v, list) or not isinstance(moduli, list):
+    if not isinstance(values, list) or not isinstance(moduli, list):
         raise ValueError("arguments to CRT_list should be lists")
-    if len(v) != len(moduli):
+    if len(values) != len(moduli):
         raise ValueError("arguments to CRT_list should be lists of the same length")
-    if not v:
+    if not values:
         return ZZ.zero()
-    if len(v) == 1:
-        return moduli[0].parent()(v[0])
-    x = v[0]
-    m = moduli[0]
+    if len(values) == 1:
+        return moduli[0].parent()(values[0])
+
+    # The result is computed using a binary tree. In typical cases,
+    # this scales much better than folding the list from one side.
     from sage.arith.functions import lcm
-    for i in range(1, len(v)):
-        x = CRT(x, v[i], m, moduli[i])
-        m = lcm(m, moduli[i])
-    return x % m
+    while len(values) > 1:
+        vs, ms = values[::2], moduli[::2]
+        for i, (v, m) in enumerate(zip(values[1::2], moduli[1::2])):
+            vs[i] = CRT(vs[i], v, ms[i], m)
+            ms[i] = lcm(ms[i], m)
+        values, moduli = vs, ms
+    return values[0] % moduli[0]
 
 
 def CRT_basis(moduli):