Chinese remainder theorem

Chinese remainder theorem
You are encouraged to solve this task according to the task description, using any language you may know.

Suppose n₁, n₂, …, n_k are positive integers that are pairwise coprime. Then, for any given sequence of integers a₁,a₂, …, a_k, there exists an integer x solving the following system of simultaneous congruences.

${\begin{aligned}x&\equiv a_{1}{\pmod {n_{1}}}\\x&\equiv a_{2}{\pmod {n_{2}}}\\&{}\ \ \vdots \\x&\equiv a_{k}{\pmod {n_{k}}}\end{aligned}}$

Furthermore, all solutions x of this system are congruent modulo the product, N = n₁n₂…n_k.

Your task is to 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 where 0 <= s < n₁n₂…n_k.

Algorithm: The following algorithm only applies if the $\scriptstyle n_{i}$'s are pairwise coprime.

Suppose, as above, that a solution is required for the system of congruences:

$x\equiv a_{i}{\pmod {n_{i}}}\quad \mathrm {for} \;i=1,\ldots ,k$

Again, to begin, the product $\scriptstyle N\;=\;n_{1}n_{2}\ldots n_{k}$ is defined. Then a solution x can be found as follows.

For each i the integers $\scriptstyle n_{i}$ and $\scriptstyle N/n_{i}$ are coprime. Using the extended Euclidean algorithm we can find integers $\scriptstyle r_{i}$ and $\scriptstyle s_{i}$ such that $\scriptstyle r_{i}n_{i}\,+\,s_{i}N/n_{i}\;=\;1$. Then, choosing the label $\scriptstyle e_{i}\;=\;s_{i}N/n_{i}$, one solution to the system of simultaneous congruences is:

$x=\sum _{i=1}^{k}a_{i}e_{i}N/n_{i}$.

Forth

Tested with GNU FORTH <lang forth>: 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 ;

</lang>

Output:

Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc.
Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license'
Type `bye' to exit
10 11  4 12  12 13  3 crt . 1000  ok
10 11  4 22   9 19  3 crt . 
:2: Invalid numeric argument
10 11  4 22   9 19  3 >>>crt<<< .

OCaml

This is using the Jane Street Ocaml Core library. <lang ocaml>open Core.Std open Option.Monad_infix

let rec egcd a b =

  if b = 0 then (1, 0)
  else
     let q = a/b and r = a mod b in
     let (s, t) = egcd b r in
        (t, s - q*t)

let mod_inv a b =

  let (x, y) = egcd a b in
     if a*x + b*y = 1 then Some x else None

let calc_inverses ns ms =

  let rec list_inverses ns ms l =
     match (ns, ms) with
        | ([], []) -> Some l
        | ([], _)
        | (_, []) -> assert false
        | (n::ns, m::ms) ->
           let inv = mod_inv n m in
              match inv with
                 | None -> None
                 | Some v -> list_inverses ns ms (v::l)
  in
     list_inverses ns ms [] >>= fun l -> Some (List.rev l)

let chinese_remainder congruences =

  let (residues, modulii) = List.unzip congruences in
  let mod_pi = List.reduce_exn modulii ~f:( * ) in
  let crt_modulii = List.map modulii ~f:(fun m -> mod_pi / m) in
  calc_inverses crt_modulii modulii >>=
     fun inverses ->
        Some (List.map3_exn residues inverses crt_modulii ~f:(fun a b c -> a*b*c)
              |> List.reduce_exn ~f:(+)
              |> fun n -> let n' = n mod mod_pi in if n' < 0 then n' + mod_pi else n')

</lang>

Output:

utop # chinese_remainder [(10, 11); (4, 12); (12, 13)];;
- : int option = Some 1000 
                                                                                                        
utop # chinese_remainder [(10, 11); (4, 22); (9, 19)];;
- : int option = None