Batched CRT, in depth: precomputed Garner for millions of values, product and remainder trees, and multi-modular algorithms

The multi-modular pattern

The pattern has four steps. First, bound the answer: |x| < B. Second, choose primes p₁..p_K that fit in a machine word, with M = product of p_i > 2B. Third, solve the whole problem modulo each prime using native integer arithmetic. Fourth, reconstruct. The third step is where the time goes, and it is independent per prime, so it spreads across cores for free. The fourth step is where careless code becomes the bottleneck or the source of wrong answers.

The factor of 2 in M > 2B matters. CRT returns x mod M, a value in [0, M). If the true answer can be negative, map values above M/2 to value minus M. That works only if every possible answer lies in (-M/2, M/2].

How many primes, and when to stop

The bound comes from the problem. For the determinant of an n by n integer matrix, Hadamard's inequality gives |det A| <= the product of the row norms. For a product of polynomials with coefficients below c in absolute value, each output coefficient is below (min(d1, d2) + 1) times c squared for degrees d1 and d2. For a dot product of length n with entries below c, the bound is n times c squared. The number of primes is K = ceil((log2 B + 1) / bits per prime). With primes just below 2³¹ (about 31 bits each) and a 10,000-bit Hadamard bound, that is ceil(10,001 / 31) = 323 primes.

Proven bounds are often far larger than the real answer. The Hadamard bound is usually very loose for random matrices. Early termination trades certainty for speed. Add primes one at a time, reconstruct incrementally, and stop when the signed result has not changed for the last few primes. If the answer were wrong, an extra random prime would have to leave it unchanged by chance, which happens with probability about 1/p. With 31-bit primes, two stable extra primes make an error extremely unlikely, but not impossible. Use the proven bound whenever the result must be certified.

Many values, few primes: Garner in batch

When K is small and fixed, as with the three or four NTT-friendly primes behind an exact convolution, the best batch method is Garner's mixed-radix form: x = d₀ + d₁m₀ + d₂m₀m₁ + ..., where each digit d_i < m_i is found using arithmetic modulo m_i only. No big integer appears until you choose to form one. Everything that depends only on the moduli is computed once: the inverses of the prefix products, and each prefix product reduced modulo every later prime. Store the residues as K rows of N values (structure of arrays), so that the inner loop runs over N values with the same constants. That loop vectorises with SIMD or NumPy and maps directly to one GPU thread per column.

from math import prod
import numpy as np

class GarnerBatch:
    # Precompute once per moduli set (each < 2**31), then reconstruct N values per call.
    def __init__(self, m):
        self.m, self.k = [int(x) for x in m], len(m)
        self.inv = [pow(prod(self.m[:i]) % self.m[i], -1, self.m[i]) for i in range(self.k)]
        self.base = [[prod(self.m[:j]) % self.m[i] for j in range(i)] for i in range(self.k)]

    def digits(self, R):                    # R: k x N uint64, row i holds x mod m[i]
        D = np.empty_like(R)
        for i, mi in enumerate(self.m):
            acc = np.zeros(R.shape[1], np.uint64)
            for j in range(i):              # value of digits 0..i-1, reduced mod m[i]
                acc = (acc + D[j] * np.uint64(self.base[i][j])) % np.uint64(mi)
            D[i] = (R[i] + np.uint64(mi) - acc) % np.uint64(mi) * np.uint64(self.inv[i]) % np.uint64(mi)
        return D                            # x = d0 + d1*m0 + d2*m0*m1 + ...

    def mod_q(self, D, q):                  # Horner in uint64; never forms x itself
        x = np.zeros(D.shape[1], np.uint64)
        for i in reversed(range(self.k)):
            x = (x * np.uint64(self.m[i] % q) + D[i] % np.uint64(q)) % np.uint64(q)
        return x

    def signed(self, D):                    # exact Python ints in (-M/2, M/2]
        M, out = prod(self.m), []
        for col in D.T:
            x = 0
            for i in reversed(range(self.k)):
                x = x * self.m[i] + int(col[i])
            out.append(x - M if x > M // 2 else x)
        return out

Check the overflow argument before you trust any reduction. Every digit and constant is below 2³¹, so each product is below 2⁶², and adding an accumulator below 2³¹ stays below 2⁶³. That is why the moduli must be under 2³¹ here. With 62-bit primes you need 128-bit products or Montgomery multiplication, as described in modular exponentiation in depth. Cost per value is about K squared over 2 multiply-reduce steps. For K = 3 that is a handful of operations per coefficient. mod_q is the step used by arbitrary-modulus convolution, covered in the number theoretic transform in depth. Note that it reduces the non-negative representative, so if the answer can be negative, subtract M mod q from values whose top digits show they exceed M/2.

Mixed-radix digits can also be compared, most significant first, with the precomputed digits of M/2, so the sign decision stays in machine words. In testing, the class above reproduced 5,004 random signed values below 2¹²⁰ across the four primes 2013265921, 1811939329, 469762049 and 998244353, both exactly and reduced modulo 10⁹+7.

One value, many primes: product and remainder trees

The other shape has one value and hundreds or thousands of moduli. Garner's O(K squared) steps now involve growing integers, and the straightforward sum of r_i times (M/m_i) times an inverse builds K numbers of the size of M. Trees fix both problems. A product tree multiplies the moduli in pairs, level by level, up to M. A remainder tree pushes x down it: each node reduces its parent's remainder by its own product, so every division is by a number of the right size. With fast multiplication, both trees cost O(M(n) log K), where M(n) is the cost of multiplying n-bit numbers. Bernstein's survey Fast Multiplication and Its Applications covers this family in depth.

For reconstruction, define c_i = r_i times the inverse of (M/m_i), all modulo m_i. The answer is the sum of c_i times M/m_i, and that sum can be built bottom-up: a node's value is its left value times the right product plus its right value times the left product. You still need (M/m_i) mod m_i for every i. Computing M mod m_i gives zero, which is useless. Instead, reduce M modulo m_i² and divide by m_i. Since M = m_i(M/m_i), the result is exactly (M/m_i) mod m_i. One remainder tree over the squares gives all K values.

def product_tree(m):
    tree = [list(m)]
    while len(tree[-1]) > 1:
        lv = tree[-1]
        tree.append([prod(lv[i:i + 2]) for i in range(0, len(lv), 2)])
    return tree                             # tree[-1][0] == prod(m)

def remainders(x, tree):
    # x mod every leaf, top-down: each node reduces its parent's remainder.
    rems = [x % tree[-1][0]]
    for lv in reversed(tree[:-1]):
        rems = [rems[i // 2] % lv[i] for i in range(len(lv))]
    return rems

def crt_tree(residues, m, tree):
    M = tree[-1][0]
    cof = [r // mi for r, mi in zip(remainders(M, product_tree([mi * mi for mi in m])), m)]
    vals = [r * pow(c, -1, mi) % mi for r, c, mi in zip(residues, cof, m)]
    for lv in tree[:-1]:                    # bottom-up: v = vL * P_R + vR * P_L
        vals = [vals[i] * lv[i + 1] + vals[i + 1] * lv[i] if i + 1 < len(lv) else vals[i]
                for i in range(0, len(lv), 2)]
    return vals[0] % M

Build the product tree once per moduli set and reuse it for every value. Measure the crossover against the naive sum on your big-integer library before committing. The same structure works for polynomials. Replace the primes with linear factors (x - a_i): the remainder tree becomes multipoint evaluation and the tree CRT becomes interpolation, which is how the fast algorithms in fast polynomial multiplication extend to evaluation.

Worked examples

Tree CRT by hand. Take moduli 3, 5, 7 and 11, so M = 1155, and residues (1, 0, 6, 10) of the unknown 1000. The product tree is [3, 5, 7, 11], [15, 77], [1155]. The cofactors M/m_i are 385, 231, 165 and 105. Reduced modulo their primes they are 1, 1, 4 and 6, with inverses 1, 1, 2 and 2. So c = (1, 0, 12 mod 7 = 5, 20 mod 11 = 9). At the first level, the node for {3, 5} is 1 times 5 + 0 times 3 = 5, and the node for {7, 11} is 5 times 11 + 9 times 7 = 118. At the root, 5 times 77 + 118 times 15 = 2155, and 2155 mod 1155 = 1000.

A multi-modular determinant. Take A = [[2, -1, 3], [4, 0, 1], [-2, 5, 6]]. The row norms are the square roots of 14, 17 and 65, so Hadamard bounds |det A| by about 124.4. We need M > 249, and the primes 17 and 19 give 323. Gaussian elimination modulo each prime gives det A as 8 mod 17 and 0 mod 19. CRT gives 76, which is below 323/2, so the answer is 76, and direct expansion confirms it. Had the residues been 9 and 0, CRT would give 247, which is above 161, and the answer would be 247 - 323 = -76. Forgetting the sign step would report 247.

Batched base conversion in RNS cryptography

Lattice cryptography keeps ciphertext coefficients in a residue number system (RNS) for its entire life. Multiplication then requires moving each coefficient from one prime basis q₁..q_k to another basis p₁..p_l. That is a batched CRT with the big integer never formed. The full-RNS variants of BFV by Bajard, Eynard, Hasan and Zucca (SAC 2016) and by Halevi, Polyakov and Shoup (CT-RSA 2019) use fast base conversion: for each target prime p_j, compute the sum over i of [x_i times (q/q_i)^-1]_qi times (q/q_i), reduced modulo p_j.

Each of the k terms is below q, so the sum equals x + αq for some integer α in [0, k): exact up to a small multiple of q. The two papers differ in how they handle that error (an extra correcting modulus versus a floating-point estimate of α). Structurally it is the Garner batch again: basis-only constants and a matrix product over N coefficients.

Failure modes

• M only just above B. Negative answers come back as large positive ones. Require M > 2B and apply the signed mapping.
• A prime that divides the problem. A determinant can be 0 modulo p while non-zero over the integers, or a needed inverse may not exist modulo p. Rational reconstruction and division-based algorithms must skip unlucky primes and use another.
• Overflow in the digit loop. 32-bit primes with uint64 accumulators can wrap silently. Prove the bound for your prime size, or use 128-bit products.
• Constants recomputed per value. Precomputation is K inversions plus O(K squared) reductions. Recomputing them for every coefficient makes the batch hundreds of times slower.
• Early termination treated as a proof. Stability over a few extra primes is probabilistic evidence. Certify with the proven bound when it matters.
• Residues out of range. Inputs that are not reduced into [0, m_i) break the subtraction step. Reduce them at the boundary.

Trade-offs

What to do next

1. Write down a proven bound B for your answer, and pick K primes with a product above 2B.
2. Precompute all modulus-only constants once, and store the residues as K rows of N values.
3. Implement GarnerBatch, then test it against Python big integers on random signed values, including M/2 and -M/2 + 1.
4. If K reaches the hundreds, switch to the product and remainder trees and reuse the tree.
5. Add an unlucky-prime check and decide whether early termination is acceptable.
6. Profile the reconstruction separately from the per-prime work. It should be a small fraction of the total.