Chinese remainder theorem

In how many ways can some of the streets around the Green Bazaar be closed so that one can still drive between any two crossings, but no loop is left? Each such set of streets is a spanning tree of the district, and in a circle with a radius of one and a half kilometres their number has 127 digits. It does not fit in a floating-point number, and computing with long numbers all the way is slow. So it is computed in parts: the remainder of division by one prime of about 67 million, by another, by a third, each time with small numbers only. The Chinese remainder theorem guarantees that from twenty-eight such remainders the number is recovered uniquely. Below this happens on real streets. Move the circle or change the radius.

Each square under the map is one “digit” of the number in a mixed radix system, given by the next modulus. While the digits are non-zero, the number is not yet assembled: there are too few moduli, and their product is smaller than the answer. Once the product of the moduli outgrows the answer, the digits become zeros: the answer is found, and further moduli only confirm it. The line “in double precision” is the same determinant computed by ordinary Gaussian elimination in floating point: of the hundred-odd digits, the first thirteen or fourteen are right.

The idea

If the moduli m₁, …, m_k are pairwise coprime, a number from 0 to m₁·…·m_k − 1 is uniquely determined by its remainders of division by each of them. So a computation that needs a large exact answer can be carried out separately modulo each small modulus, with machine numbers, without long arithmetic, even on different computers, and at the end the answer is assembled from the residues. Garner's algorithm assembles it: the answer is written as c₀ + c₁·m₁ + c₂·m₁m₂ + …, and each “digit” c_i is found modulo m_i+1 in one pass over those already found.

The theorem itself and its classical proof by a table of remainders are taken apart in the chapter “Sunzi’s problem” of “Mathematics, the Queen of the Sciences”. The same place has the formula that assembles the answer from “parts”: numbers equal to one modulo one modulus and zero modulo the others. Here we are interested in the theorem as a tool: how to recover a number quickly, how many moduli are needed and where it pays off.

Garner's algorithm

The formula with parts requires multiplying large numbers: the part for the modulus m_i is the product of all the other moduli. Harvey Garner in 1959, working on arithmetic in residue number systems for computers, proposed assembling the number differently, by the digits of a mixed radix system. Write

Modulo m₁ all the terms except the first vanish: c₀ = r₁. Modulo m₂ all the terms from the third on vanish, so r₂ ≡ c₀ + c₁m₁, and c₁ = (r₂ − c₀)·m₁⁻¹ mod m₂. And so on: each next digit is computed modulo its own modulus from those already known, and the inverses come from the extended Euclidean algorithm. All numbers in these computations are smaller than the largest modulus; long arithmetic is needed only at the very end, to add up the c_i multiplied by the products of the moduli.

See what this looks like on a calendar. A year of the Julian period, a number from 1 to 7980, is given by three remainders: the indiction modulo 15, the golden number modulo 19 and the solar cycle modulo 28. Each row of the scene is one of Garner's digits: jumps by 1, then by 15, then by 15 · 19 = 285. A jump by the product of the previous moduli does not change the previous remainders, and that is the whole algorithm.

The digits c_i have a handy property, visible in the first scene: if the answer is smaller than m₁·…·m_j, all digits from c_j on are zeros. So the computation can be stopped when several new moduli in a row have given zero digits. Strictly this proves nothing, since a zero digit turns up by chance with probability 1/m, and if a guarantee is needed, the number of moduli is chosen from a bound on the answer known in advance.

Trees, the determinant and how many moduli are needed

Where the determinant in the first scene comes from. In 1847 Gustav Kirchhoff, calculating currents in electrical circuits, found a formula for the number of spanning trees of a graph. Take the matrix with the degrees of the vertices on the diagonal and, for each edge uv, minus one in the cells (u, v) and (v, u), and cross out one row and the same column. The determinant of what is left equals the number of spanning trees. The proof of this formula is longer than the article allows; it can be found in graph theory textbooks under the name “matrix tree theorem”. The same matrix, with conductances instead of ones, is the system of equations from the article on Gaussian elimination, where current flows through the streets of Almaty.

The determinant of an integer matrix is an integer, and modulo a prime p it is computed by the same Gaussian elimination, with division replaced by multiplication by the inverse modulo p. To recover the answer, the product of the moduli must exceed the answer itself, so a bound known in advance is needed. For trees it comes out very simply.

The bound is crude. For the district by the bazaar it promises up to 213 digits, while the answer has 127; in the scene this shows as the tail of zero Garner digits. But it is honest: how many moduli to take is known before the computation starts. For arbitrary matrices Hadamard's inequality plays the same role: the determinant is at most the product of the lengths of the rows.

One more subtlety is an “unlucky” modulus. Gaussian elimination in the scene does not swap rows: Kirchhoff's matrix with a row and column crossed out is positive definite, and in ordinary numbers a pivot is never zero. Modulo p it can be, with probability about 1/p per step. Such a modulus is simply skipped and the next one taken: the residues from the other primes are enough.

Code

The program takes from the map a district with a radius of one kilometre around the Green Bazaar, numbers the crossings so that neighbours get close numbers, so that the non-zero entries of the matrix lie in a narrow band along the diagonal, and computes the determinant modulo primes below 2³¹. The product of two residues is below 2⁶², so in Go uint64 is enough, and math/big is needed only at the very end, in Garner's algorithm. The last line is a check by a modulus that took no part in the computation.

district: 237 crossings, 328 streets, matrix 236 x 236, band 38
bound: at most 97 digits
11 primes below 2^31, the largest 2147483647, the smallest 2147483423
spanning trees: 34613761769290138902561752769643080652778611353227531802400
59 digits
check modulo 1000000007: 525194220 = 525194220

The second program is the calendar and the general case. The Julian period is recovered in two ways: by Garner's algorithm and by the formula that Jacques de Billy published in 1665 and Gauss gave in 1801: the year of the period equals (6916·i + 4200·g + 4845·s) mod 7980. The numbers 6916, 4200 and 4845 are those very “parts”: 6916 = 19 · 28 · 13 leaves remainder 1 on division by 15 and is divisible by 19 and 28, and so on.

year       indiction golden solar   Julian period  de Billy  mixed radix
2026 AD            4     13    19            6739      6739  4 + 12*15 + 23*285
1583 AD           11      7    24            6296      6296  11 + 1*15 + 22*285
1 AD               4      2    10            4714      4714  4 + 10*15 + 16*285
1 BC               3      1     9            4713      4713  3 + 10*15 + 16*285
4713 BC            1      1     1               1         1  1 + 0*15 + 0*285

If the moduli are not coprime, a solution does not always exist. The congruences x ≡ r₁ (mod m₁) and x ≡ r₂ (mod m₂) are compatible exactly when r₁ ≡ r₂ modulo g = gcd(m₁, m₂), and then the solution is unique modulo lcm(m₁, m₂). The function merge finds it with the extended Euclidean algorithm: x = r₁ + m₁t, where t solves (m₁/g)·t ≡ (r₂ − r₁)/g modulo m₂/g. The two congruences x ≡ 2 (mod 4) and x ≡ 3 (mod 6) are incompatible, since the first requires an even x and the second an odd one, while x ≡ 2 (mod 4) and x ≡ 4 (mod 6) give x ≡ 10 (mod 12).

What it costs

Recovery from k moduli by Garner's algorithm is O(k²) operations on machine numbers plus O(k²) for the final assembly of the long number. That is little: the main work is the computations modulo each modulus, k of them, and they do not depend on each other, so they can be handed out to different cores or machines. The number of moduli grows linearly with the length of the answer: one 31-bit modulus for every 9 decimal digits.

Does it pay off compared with exact computation in fractions? A measurement: the same determinant for districts of different sizes, computed two ways, by moduli and by ordinary Gaussian elimination on fractions that reduce themselves.

In Go the moduli win 9–15 times on matrices from seventy rows: on a matrix of 457 rows, 62 milliseconds against 586. In Python it is almost a draw, and on large matrices fractions are even slightly faster. The reason is the language, not the algorithm: the long-number arithmetic inside Fraction is written in C, while the loop over moduli is in Python, and its overhead eats up the gain from small numbers. It is the same story as with built-in multiplication in the article on Karatsuba multiplication: counting operations is useful, but the measurement decides.

Where it is used

• Exact linear algebra. The FLINT library computes the determinant of an integer matrix from 25 rows over many moduli, and smaller matrices by Bareiss's method, without fractions. The same technique works for polynomials: the greatest common divisor and factorisation are easier to find modulo primes and then lift the result.
• RSA. The private operation c^d mod N is computed separately modulo p and q, with numbers and exponents half as long, and glued together by the Chinese theorem. Exponentiation costs roughly the cube of the length, each half an eighth of the whole, and both together about a quarter; that is why an RSA private key stores not only d but also p, q, d mod (p − 1), d mod (q − 1) and q⁻¹ mod p. The speed-up has a price: if a hardware fault spoils one of the halves, the wrong signature gives away a factor of N through a gcd, so implementations check a signature before handing it out.
• Multiplying long numbers. The number-theoretic transform multiplies polynomials exactly, but modulo a prime. If the coefficients of the product are larger than the modulus, three primes of the form c·2^k + 1 are taken, for example 998,244,353 = 119·2²³ + 1, and the three results are glued together by Garner's algorithm.
• Residue number systems. In a representation by residues, addition and multiplication go independently modulo each modulus, without carries between digits. This is what the paper by Garner was written for: such systems were studied for computers and are used in specialised digital signal processing hardware. Comparing numbers and dividing in residues is hard, though, and that is what conversion to ordinary notation is needed for.
• Chronology. The Julian period of the second scene was devised by Joseph Scaliger in 1583 to bring the dates of different calendars to one count; its first year, 4713 BC, is the one where all three cycles equal one. Since 1849, after John Herschel's book, astronomers have counted Julian days from noon on 1 January of that year, and they still do.

Where it trips up

• Moduli with a common divisor. The formula with parts and Garner's algorithm require coprime moduli: 4 has no inverse modulo 6. The general case needs a compatibility check and gluing through the lcm, as in the function merge.
• Too few moduli. If the product of the moduli is smaller than the answer, Garner will honestly return the remainder of the answer modulo this product: a plausible and wrong number. The bound on the answer is needed before the computation.
• A negative answer. A determinant can be negative, while residues cannot. The recovered x from the range [0, M) must be read as a number from [−M/2, M/2): if x > M/2, the answer is x − M. One more modulus is then needed.
• Overflow. In Go and C the product of two residues modulo a number around 2⁶³ does not fit in 64 bits. Moduli up to 2³¹ are the simplest way to do without 128-bit multiplication.

History

The problem of objects that, counted by threes, fives and sevens, leave remainders 2, 3 and 2 is in the treatise Sunzi Suanjing, written in the 3rd–5th centuries; a general method for solving such systems was set out by Qin Jiushao in 1247. In the language of congruences the theorem was stated and proved by Gauss in the Disquisitiones Arithmeticae of 1801, where, among the examples, is the Julian period. As an algorithm for computers it came back in 1959 in Harvey Garner's paper on residue number systems. Later, multi-modular computation became a standard technique of computer algebra: this is how the greatest common divisor of polynomials is found, where exact coefficients swell, and how determinants are computed.

Problems

The first is a calendar one: recover the year from three cycles without being off by one at the boundary of eras. In the second you count the spanning trees of a graph exactly, when their number has long stopped fitting in any machine type; one of the tests is a real district of Almaty.