compute list CRT via tree

[yyyyx4]

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: 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)

Component: algebra

Author: Lorenz Panny

Branch/Commit: 5f3a23f

Reviewer: Kwankyu Lee

Issue created by migration from https://trac.sagemath.org/ticket/34512

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

eac0595

avoid mutating inputs by doing first layer out-of-place

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 63f6f3c to eac0595

[dimpase]

comment:5

test, please ignore

[kwankyu]

comment:6

I simplified(?) the code

--- a/src/sage/arith/misc.py
+++ b/src/sage/arith/misc.py
@@ -3302,13 +3302,14 @@ def CRT_list(v, moduli):
     if len(v) == 1:
         return moduli[0].parent()(v[0])
     from sage.arith.functions import lcm
-    v = [CRT(v[i], v[i+1], moduli[i], moduli[i+1]) for i in range(0, len(v)-1, 2)] + v[len(v)//2*2:]
-    moduli = [lcm(moduli[i], moduli[i+1]) for i in range(0, len(moduli)-1, 2)] + moduli[len(moduli)//2*2:]
     while len(v) > 1:
-        for i in range(0, len(v)-1, 2):
-            v[i] = CRT(v[i], v[i+1], moduli[i], moduli[i+1])
-            moduli[i] = lcm(moduli[i], moduli[i+1])
-        v, moduli = v[::2], moduli[::2]
+        vs = [CRT(v[i], v[i + 1], moduli[i], moduli[i + 1]) for i in range(0, len(v) - 1, 2)]
+        ms = [lcm(moduli[i], moduli[i + 1]) for i in range(0, len(v) - 1, 2)]
+        if len(v) % 2:
+            vs.append(v[-1])
+            ms.append(moduli[-1])
+        v = vs
+        moduli = ms
     return v[0] % moduli[0]

This is slightly slower than your code but easier to read. What do you think?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

65bce21

improve readability

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from eac0595 to 65bce21

[yyyyx4]

comment:8

Thanks! I continued tweaking it and ended up with this. It seems to be just as fast as my first version.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 65bce21 to 819d3d3

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

ef41957	Merge branch 'develop' into t/34512/public/compute_list_CRT_via_tree
819d3d3	Move implementation comment to code block

[kwankyu]

comment:10

Thanks. Looks nice.

I moved the implementation comment to the code block. Other than that, I am positive.

[kwankyu]

Reviewer: Kwankyu Lee

[vbraun]

comment:11

sage -t --warn-long 46.9 --random-seed=47449161950200526044278063017024899534 src/sage/rings/polynomial/polynomial_quotient_ring.py
**********************************************************************
File "src/sage/rings/polynomial/polynomial_quotient_ring.py", line 1643, in sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.S_units
Failed example:
    [u for u, o in L.S_units([]) if o is Infinity]
Exception raised:
    Traceback (most recent call last):
      File "/home/release/Sage/src/sage/doctest/forker.py", line 695, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/release/Sage/src/sage/doctest/forker.py", line 1093, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.rings.polynomial.polynomial_quotient_ring.PolynomialQuotientRing_generic.S_units[11]>", line 1, in <module>
        [u for u, o in L.S_units([]) if o is Infinity]
      File "/home/release/Sage/src/sage/rings/polynomial/polynomial_quotient_ring.py", line 1687, in S_units
        poly_unit = self(sage.arith.all.crt(prod_unit, moduli))
      File "/home/release/Sage/src/sage/arith/misc.py", line 3203, in crt
        return CRT_list(a, b)
      File "/home/release/Sage/src/sage/arith/misc.py", line 3320, in CRT_list
        return values[0] % moduli[0]
      File "sage/rings/polynomial/polynomial_element.pyx", line 2882, in sage.rings.polynomial.polynomial_element.Polynomial.__mod__
        _, R = self.quo_rem(other)
      File "sage/structure/element.pyx", line 4500, in sage.structure.element.coerce_binop.new_method
        a, b = coercion_model.canonical_coercion(self, other)
      File "sage/structure/coerce.pyx", line 1393, in sage.structure.coerce.CoercionModel.canonical_coercion
        raise TypeError("no common canonical parent for objects with parents: '%s' and '%s'"%(xp, yp))
    TypeError: no common canonical parent for objects with parents: 'Univariate Polynomial Ring in x over Number Field in a with defining polynomial x^2 + 3 with a = 1.732050807568878?*I' and 'Univariate Polynomial Ring in y over Number Field in a with defining polynomial x^2 + 3 with a = 1.732050807568878?*I'

and others

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 819d3d3 to 5f3a23f

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

432493f	Merge branch 'develop' into t/34512/public/compute_list_CRT_via_tree
5f3a23f	Recover the case for single pair

[kwankyu]

comment:14

Passed BUILD&TEST.

[vbraun]

Changed branch from public/compute_list_CRT_via_tree to 5f3a23f