Problem Statement

Suppose

n

1

{\displaystyle n_{1}}

,

n

2

{\displaystyle n_{2}}

,

…

{\displaystyle \ldots }

,

n

k

{\displaystyle n_{k}}

are positive integers that are pairwise co-prime.   Then, for any given sequence of integers

a

1

{\displaystyle a_{1}}

,

a

2

{\displaystyle a_{2}}

,

…

{\displaystyle \dots }

,

a

k

{\displaystyle a_{k}}

,   there exists an integer

x

{\displaystyle x}

solving the following system of simultaneous congruences: Furthermore, all solutions

x

{\displaystyle x}

of this system are congruent modulo the product,

N

n

1

n

2

…

n

k

{\displaystyle N=n_{1}n_{2}\ldots n_{k}}

.

Write a program to solve a system of linear congruences by applying the   Chinese Remainder Theorem. If the system of equations cannot be solved, your program must somehow indicate this. (It may throw an exception or return a special false value.) Since there are infinitely many solutions, the program should return the unique solution

s

{\displaystyle s}

where

0 ≤ s ≤

n

1

n

2

…

n

k

{\displaystyle 0\leq s\leq n_{1}n_{2}\ldots n_{k}}

.

Show the functionality of this program by printing the result such that the

n

{\displaystyle n}

's   are

[ 3 , 5 , 7 ]

{\displaystyle [3,5,7]}

and the

a

{\displaystyle a}

's   are

[ 2 , 3 , 2 ]

{\displaystyle [2,3,2]}

.

Algorithm:   The following algorithm only applies if the

n

i

{\displaystyle n_{i}}

's   are pairwise co-prime. Suppose, as above, that a solution is required for the system of congruences: Again, to begin, the product

N

n

1

n

2

…

n

k

{\displaystyle N=n_{1}n_{2}\ldots n_{k}}

is defined. Then a solution

x

{\displaystyle x}

can be found as follows: For each

i

{\displaystyle i}

,   the integers

n

i

{\displaystyle n_{i}}

and

N

/

n

i

{\displaystyle N/n_{i}}

are co-prime. Using the   Extended Euclidean algorithm,   we can find integers

r

i

{\displaystyle r_{i}}

and

s

i

{\displaystyle s_{i}}

such that

r

i

n

i

s

i

N

/

n

i

= 1

{\displaystyle r_{i}n_{i}+s_{i}N/n_{i}=1}

. Then, one solution to the system of simultaneous congruences is: and the minimal solution,

Let's start with the solution:

: egcd ( a b -- a b )
    dup 0= IF
        2drop 1 0
    ELSE
        dup -rot /mod               \ -- b r=a%b q=a/b
        -rot recurse                \ -- q (s,t) = egcd(b, r)
        >r swap r@ * - r> swap      \ -- t (s - q*t)
    THEN ;

: egcd>gcd  ( a b x y -- n )  \ calculate gcd from egcd
    rot * -rot * + ;

: mod-inv  ( a m -- a' )     \ modular inverse with coprime check
    2dup egcd over >r egcd>gcd r> swap 1 <> -24 and throw ;

: array-product ( adr count -- n )
    1 -rot  cells bounds ?DO  i @ *  cell +LOOP ;

: crt-from-array  ( adr1 adr2 count -- n )
    2dup array-product   locals| M count m[] a[] |
    0  \ result
    count 0 DO
        m[] i cells + @
        dup M swap /
        dup rot mod-inv *
        a[] i cells + @ * +
    LOOP  M mod ;

create crt-residues[]  10 cells allot
create crt-moduli[]    10 cells allot

: crt ( .... n -- n )  \ takes pairs of "n (mod m)" from stack.
    10 min  locals| n |
    n 0 DO
        crt-moduli[] i cells + !
        crt-residues[] i cells + !
    LOOP
    crt-residues[] crt-moduli[] n crt-from-array ;