/** 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
}