Chinese remainder theorem
You are encouraged to solve this task according to the task description, using any language you may know.
Suppose , , , are positive integers that are pairwise coprime. Then, for any given sequence of integers , , , , there exists an integer solving the following system of simultaneous congruences.
Furthermore, all solutions of this system are congruent modulo the product, .
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 where .
Show the functionality of this program by printing the result such that the 's are and the 's are .
Algorithm: The following algorithm only applies if the 's are pairwise coprime.
Suppose, as above, that a solution is required for the system of congruences:
Again, to begin, the product is defined. Then a solution can be found as follows.
For each , the integers and are coprime. Using the extended Euclidean algorithm we can find integers and such that . Then, choosing the label , one solution to the system of simultaneous congruences is:
- .
C
When n are not pariwise coprime, the program crashes due to division by zero, which is one way to convey error. <lang c>#include <stdio.h>
// returns x where (a * x) % b == 1 int mul_inv(int a, int b) { int b0 = b, t, q; int x0 = 0, x1 = 1; if (b == 1) return 1; while (a > 1) { q = a / b; t = b, b = a % b, a = t; t = x0, x0 = x1 - q * x0, x1 = t; } if (x1 < 0) x1 += b0; return x1; }
int chinese_remainder(int *n, int *a, int len) { int p, i, prod = 1, sum = 0;
for (i = 0; i < len; i++) prod *= n[i];
for (i = 0; i < len; i++) { p = prod / n[i]; sum += a[i] * mul_inv(p, n[i]) * p; }
return sum % prod; }
int main(void) { int n[] = { 3, 5, 7 }; int a[] = { 2, 3, 2 };
printf("%d\n", chinese_remainder(n, a, sizeof(n)/sizeof(n[0]))); return 0; }</lang>
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
Perl 6
<lang perl6># returns x where (a * x) % b == 1 sub mul-inv($a is copy, $b is copy) {
return 1 if $b == 1; my ($b0, @x) = $b, 0, 1; ($a, $b, @x) = (
$b, $a % $b, @x[1] - ($a div $b)*@x[0], @x[0]
) while $a > 1; @x[1] += $b0 if @x[1] < 0; return @x[1];
}
sub chinese-remainder(*@n) {
my \N = [*] @n; -> *@a {
N R% [+] map { my \p = N div @n[$_]; @a[$_] * mul-inv(p, @n[$_]) * p }, ^@n
}
}
say chinese-remainder(3, 5, 7)(2, 3, 2);</lang>
- Output:
23
Tcl
<lang tcl>proc ::tcl::mathfunc::mulinv {a b} {
if {$b == 1} {return 1} set b0 $b; set x0 0; set x1 1 while {$a > 1} {
set x0 [expr {$x1 - ($a / $b) * [set x1 $x0]}] set b [expr {$a % [set a $b]}]
} incr x1 [expr {($x1 < 0) * $b0}]
} proc chineseRemainder {nList aList} {
set sum 0; set prod [::tcl::mathop::* {*}$nList] foreach n $nList a $aList {
set p [expr {$prod / $n}] incr sum [expr {$a * mulinv($p, $n) * $p}]
} expr {$sum % $prod}
} puts [chineseRemainder {3 5 7} {2 3 2}]</lang>
- Output:
23