/** This contains a function to solve a system of modular integer congruences using the Chinese Remainder Theorem. See: https://brilliant.org/wiki/chinese-remainder-theorem/ Given pairwise coprime positive integer moduli m1, m2, m3, ... mk and arbitrary integer remainders r1, r2, r3, ... rk, this solves the system of simultaneous congruences x ≡ r1 (mod m1) x ≡ r2 (mod m2) x ≡ r3 (mod m3) ... x ≡ rk (mod mk) And returns the unique solution (x mod N) where N is the product of all the m terms. arguments: [r, m] where r and m are arrays of the remainder terms r and the modulus terms m respectively. These must be of the same length. returns x, the unique solution mod N where N is the product of all the M terms. Example: ChineseRemainder[[1,4,6], [3,5,7]] == 34 The "comet" puzzle on the Brilliant site above can be solved by: ChineseRemainder[[2017,2014,2008],[3,8,13],2017] == 2086 Note: This does not verify that the moduli are co-prime. But maybe that can be done? See https://en.wikipedia.org/wiki/Chinese_remainder_theorem#Generalization_to_non-coprime_moduli This does not currently seem to work for negative numbers. Fix that. */ ChineseRemainder[r, m] := { return ChineseRemainder[r, m, 0] } /** Returns a solution using the Chinese Remainder Theorem as above, but making the solution x be the smallest solution >= d. */ ChineseRemainder[r, m, d] := { if length[r] != length[m] { println["ChineseRemainder: r and m must be arrays of the same length."] return undef } N = product[m] y = new array z = new array x = 0 for i = rangeOf[m] { y@i = N / m@i z@i = modInverse[y@i, m@i] if z@i == undef { println["ChineseRemainder: modInverse returned undef for modInverse[" + y@i + ", " + m@i + "]"] return undef } x = x + r@i y@i z@i } xp = x mod N f = d div N r = f * N + xp if r < d r = r + N return r }