Chinese remainder theorem

From Rosetta Code
Task
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 co-prime.  

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,   .


Task

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 co-prime.

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 co-prime.

Using the   Extended Euclidean algorithm,   we can find integers     and     such that   .

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

and the minimal solution,

.



360 Assembly[edit]

Translation of: REXX
*        Chinese remainder theorem 06/09/2015
CHINESE CSECT
USING CHINESE,R12 base addr
LR R12,R15
BEGIN LA R9,1 m=1
LA R6,1 j=1
LOOPJ C R6,NN do j=1 to nn
BH ELOOPJ
LR R1,R6 j
SLA R1,2 j*4
M R8,N-4(R1) m=m*n(j)
LA R6,1(R6) j=j+1
B LOOPJ
ELOOPJ LA R6,1 x=1
LOOPX CR R6,R9 do x=1 to m
BH ELOOPX
LA R7,1 i=1
LOOPI C R7,NN do i=1 to nn
BH ELOOPI
LR R1,R7 i
SLA R1,2 i*4
LR R5,R6 x
LA R4,0
D R4,N-4(R1) x//n(i)
C R4,A-4(R1) if x//n(i)^=a(i)
BNE ITERX then iterate x
LA R7,1(R7) i=i+1
B LOOPI
ELOOPI MVC PG(2),=C'x='
XDECO R6,PG+2 edit x
XPRNT PG,14 print buffer
B RETURN
ITERX LA R6,1(R6) x=x+1
B LOOPX
ELOOPX XPRNT NOSOL,17 print
RETURN XR R15,R15 rc=0
BR R14
NN DC F'3'
N DC F'3',F'5',F'7'
A DC F'2',F'3',F'2'
PG DS CL80
NOSOL DC CL17'no solution found'
YREGS
END CHINESE
Output:
x=          23

Ada[edit]

Using the package Mod_Inv from [[1]].

with Ada.Text_IO, Mod_Inv;
 
procedure Chin_Rema is
N: array(Positive range <>) of Positive := (3, 5, 7);
A: array(Positive range <>) of Positive := (2, 3, 2);
Tmp: Positive;
Prod: Positive := 1;
Sum: Natural := 0;
 
begin
for I in N'Range loop
Prod := Prod * N(I);
end loop;
 
for I in A'Range loop
Tmp := Prod / N(I);
Sum := Sum + A(I) * Mod_Inv.Inverse(Tmp, N(I)) * Tmp;
end loop;
Ada.Text_IO.Put_Line(Integer'Image(Sum mod Prod));
end Chin_Rema;

Bracmat[edit]

Translation of: C
( ( mul-inv
= a b b0 q x0 x1
.  !arg:(?a.?b:?b0)
& ( !b:1
| 0:?x0
& 1:?x1
& whl
' ( !a:>1
& (!b.mod$(!a.!b):?q.!x1+-1*!q*!x0.!x0)
 : (?a.?b.?x0.?x1)
)
& ( !x1:<0&!b0+!x1
| !x1
)
)
)
& ( chinese-remainder
= n a as p ns ni prod sum
.  !arg:(?n.?a)
& 1:?prod
& 0:?sum
& !n:?ns
& whl'(!ns:%?ni ?ns&!prod*!ni:?prod)
& !n:?ns
& !a:?as
& whl
' ( !ns:%?ni ?ns
& !as:%?ai ?as
& div$(!prod.!ni):?p
& !sum+!ai*mul-inv$(!p.!ni)*!p:?sum
)
& mod$(!sum.!prod):?arg
& !arg
)
& 3 5 7:?n
& 2 3 2:?a
& put$(str$(chinese-remainder$(!n.!a) \n))
);

Output:

23

C[edit]

When n are not pairwise coprime, the program crashes due to division by zero, which is one way to convey error.

#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;
}

Common Lisp[edit]

Using function invmod from [[2]].

 
(defun chinese-remainder (am)
"Calculates the Chinese Remainder for the given set of integer modulo pairs.
Note: All the ni and the N must be coprimes."

(loop :for (a . m) :in am
:with mtot = (reduce #'* (mapcar #'(lambda(X) (cdr X)) am))
:with sum = 0
:finally (return (mod sum mtot))
:do
(incf sum (* a (invmod (/ mtot m) m) (/ mtot m)))))
 

Output:

* (chinese-remainder '((2 . 3) (3 . 5) (2 . 7)))

23
* (chinese-remainder '((10 . 11) (4 . 12) (12 . 13)))

1000
* (chinese-remainder '((19 . 100) (0 . 23)))

1219
* (chinese-remainder '((10 . 11) (4 . 22) (9 . 19)))

debugger invoked on a SIMPLE-ERROR in thread
#<THREAD "main thread" RUNNING {1002A8B1B3}>:
  invmod: Values 418 and 11 are not coprimes.

Type HELP for debugger help, or (SB-EXT:EXIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
  0: [ABORT] Exit debugger, returning to top level.

(INVMOD 418 11)
0] 

D[edit]

Translation of: Python
import std.stdio, std.algorithm;
 
T chineseRemainder(T)(in T[] n, in T[] a) pure nothrow @safe @nogc
in {
assert(n.length == a.length);
} body {
static T mulInv(T)(T a, T b) pure nothrow @safe @nogc {
auto b0 = b;
T x0 = 0, x1 = 1;
if (b == 1)
return T(1);
while (a > 1) {
immutable q = a / b;
immutable amb = a % b;
a = b;
b = amb;
immutable xqx = x1 - q * x0;
x1 = x0;
x0 = xqx;
}
if (x1 < 0)
x1 += b0;
return x1;
}
 
immutable prod = reduce!q{a * b}(T(1), n);
 
T p = 1, sm = 0;
foreach (immutable i, immutable ni; n) {
p = prod / ni;
sm += a[i] * mulInv(p, ni) * p;
}
return sm % prod;
}
 
void main() {
immutable n = [3, 5, 7],
a = [2, 3, 2];
chineseRemainder(n, a).writeln;
}
Output:
23

EchoLisp[edit]

egcd - extended gcd - and crt-solve - chinese remainder theorem solve - are included in math.lib.

 
(lib 'math)
math.lib v1.10 ® EchoLisp
Lib: math.lib loaded.
 
(crt-solve '(2 3 2) '(3 5 7))
23
(crt-solve '(2 3 2) '(7 1005 15))
💥 error: mod[i] must be co-primes : assertion failed : 1005
 

Elixir[edit]

Translation of: Ruby

Brute-force:

defmodule Chinese do
def remainder(mods, remainders) do
max = Enum.reduce(mods, fn x,acc -> x*acc end)
Enum.zip(mods, remainders)
|> Enum.map(fn {m,r} -> Enum.take_every(r..max, m) |> Enum.into(HashSet.new) end)
|> Enum.reduce(fn set,acc -> Set.intersection(set, acc) end)
|> Set.to_list
end
end
 
IO.inspect Chinese.remainder([3,5,7], [2,3,2])
IO.inspect Chinese.remainder([10,4,9], [11,22,19])
IO.inspect Chinese.remainder([11,12,13], [10,4,12])
Output:
[23]
[]
[1000]

Erlang[edit]

Translation of: OCaml
-module(crt).
-import(lists, [zip/2, unzip/1, foldl/3, sum/1]).
-export([egcd/2, mod/2, mod_inv/2, chinese_remainder/1]).
 
egcd(_, 0) -> {1, 0};
egcd(A, B) ->
{S, T} = egcd(B, A rem B),
{T, S - (A div B)*T}.
 
mod_inv(A, B) ->
{X, Y} = egcd(A, B),
if
A*X + B*Y =:= 1 -> X;
true -> undefined
end.
 
mod(A, M) ->
X = A rem M,
if
X < 0 -> X + M;
true -> X
end.
 
calc_inverses([], []) -> [];
calc_inverses([N | Ns], [M | Ms]) ->
case mod_inv(N, M) of
undefined -> undefined;
Inv -> [Inv | calc_inverses(Ns, Ms)]
end.
 
chinese_remainder(Congruences) ->
{Residues, Modulii} = unzip(Congruences),
ModPI = foldl(fun(A, B) -> A*B end, 1, Modulii),
CRT_Modulii = [ModPI div M || M <- Modulii],
case calc_inverses(CRT_Modulii, Modulii) of
undefined -> undefined;
Inverses ->
Solution = sum([A*B || {A,B} <- zip(CRT_Modulii,
[A*B || {A,B} <- zip(Residues, Inverses)])]),
mod(Solution, ModPI)
end.
Output:
16> crt:chinese_remainder([{10,11}, {4,12}, {12,13}]).
1000
17> crt:chinese_remainder([{10,11}, {4,22}, {9,19}]).
undefined
18> crt:chinese_remainder([{2,3}, {3,5}, {2,7}]).
23

Forth[edit]

Tested with GNU 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 ;
 
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<<< .

FunL[edit]

import integers.modinv
 
def crt( congruences ) =
N = product( n | (_, n) <- congruences )
sum( a*modinv(N/n, n)*N/n | (a, n) <- congruences ) mod N
 
println( crt([(2, 3), (3, 5), (2, 7)]) )
Output:
23

Go[edit]

Go has the Extended Euclidean algorithm in the GCD function for big integers in the standard library. GCD will return 1 only if numbers are coprime, so a result != 1 indicates the error condition.

package main
 
import (
"fmt"
"math/big"
)
 
var one = big.NewInt(1)
 
func crt(a, n []*big.Int) (*big.Int, error) {
p := new(big.Int).Set(n[0])
for _, n1 := range n[1:] {
p.Mul(p, n1)
}
var x, q, s, z big.Int
for i, n1 := range n {
q.Div(p, n1)
z.GCD(nil, &s, n1, &q)
if z.Cmp(one) != 0 {
return nil, fmt.Errorf("%d not coprime", n1)
}
x.Add(&x, s.Mul(a[i], s.Mul(&s, &q)))
}
return x.Mod(&x, p), nil
}
 
func main() {
n := []*big.Int{
big.NewInt(3),
big.NewInt(5),
big.NewInt(7),
}
a := []*big.Int{
big.NewInt(2),
big.NewInt(3),
big.NewInt(2),
}
fmt.Println(crt(a, n))
}
Output:

Two values, the solution x and an error value.

23 <nil>

Haskell[edit]

Translation of: Erlang
module CRT ( chineseRemainder ) where
 
egcd :: Integral a => a -> a -> (a, a)
egcd _ 0 = (1, 0)
egcd a b = (t, s - q * t)
where
(s, t) = egcd b r
(q, r) = a `quotRem` b
 
modInv :: Integral a => a -> a -> Maybe a
modInv a b = case egcd a b of
(x, y) | a * x + b * y == 1 -> Just x
| otherwise -> Nothing
 
chineseRemainder :: Integral a => [a] -> [a] -> Maybe a
chineseRemainder residues modulii = do
inverses <- sequence $ zipWith modInv crtModulii modulii
return . (`mod` modPI) . sum $
zipWith (*) crtModulii (zipWith (*) residues inverses)
where
modPI = product modulii
crtModulii = map (modPI `div`) modulii
 
Output:
λ> chineseRemainder [10, 4, 12] [11, 12, 13]
Just 1000
λ> chineseRemainder [10, 4, 9] [11, 22, 19]
Nothing
λ> chineseRemainder [2, 3, 2] [3, 5, 7]
Just 23

Icon and Unicon[edit]

Translation of: Python
with error check added.

Works in both languages:

link numbers   # for gcd()
 
procedure main()
write(cr([3,5,7],[2,3,2]) | "No solution!")
write(cr([10,4,9],[11,22,19]) | "No solution!")
end
 
procedure cr(n,a)
if 1 ~= gcd(n[i := !*n],a[i]) then fail # Not pairwise coprime
(prod := 1, sm := 0)
every prod *:= !n
every p := prod/(ni := n[i := !*n]) do sm +:= a[i] * mul_inv(p,ni) * p
return sm%prod
end
 
procedure mul_inv(a,b)
if b = 1 then return 1
(b0 := b, x0 := 0, x1 := 1)
while q := (1 < a)/b do {
(t := a, a := b, b := t%b)
(t := x0, x0 := x1-q*t, x1 := t)
}
return if x1 < 0 then x1+b0 else x1
end

Output:

->crt
23
No solution!
->

J[edit]

Solution (brute force):
   crt =: (1 + ] - {:@:[ -: {.@:[ | ])^:_&0@:,:
Example:
   3 5 7 crt 2 3 2 
23
11 12 13 crt 10 4 12
1000

Notes: This is a brute force approach and does not meet the requirement for explicit notification of an an unsolvable set of equations (it just spins forever). A much more thorough and educational approach can be found on the J wiki's Essay on the Chinese Remainder Thereom.

Java[edit]

Translation of Python via D

Works with: Java version 8
import static java.util.Arrays.stream;
 
public class ChineseRemainderTheorem {
 
public static int chineseRemainder(int[] n, int[] a) {
 
int prod = stream(n).reduce(1, (i, j) -> i * j);
 
int p, sm = 0;
for (int i = 0; i < n.length; i++) {
p = prod / n[i];
sm += a[i] * mulInv(p, n[i]) * p;
}
return sm % prod;
}
 
private static int mulInv(int a, int b) {
int b0 = b;
int x0 = 0;
int x1 = 1;
 
if (b == 1)
return 1;
 
while (a > 1) {
int q = a / b;
int amb = a % b;
a = b;
b = amb;
int xqx = x1 - q * x0;
x1 = x0;
x0 = xqx;
}
 
if (x1 < 0)
x1 += b0;
 
return x1;
}
 
public static void main(String[] args) {
int[] n = {3, 5, 7};
int[] a = {2, 3, 2};
System.out.println(chineseRemainder(n, a));
}
}
23

jq[edit]

This implementation is similar to the one in C, but raises an error if there is no solution, as illustrated in the last example.

# mul_inv(a;b) returns x where (a * x) % b == 1, or else null
def mul_inv(a; b):
 
# state: [a, b, x0, x1]
def iterate:
.[0] as $a | .[1] as $b
| if $a > 1 then
if $b == 0 then null
else ($a / $b | floor) as $q
| [$b, ($a % $b), (.[3] - ($q * .[2])), .[2]] | iterate
end
else .
end ;
 
if (b == 1) then 1
else [a,b,0,1] | iterate
| if . == null then .
else .[3] | if . < 0 then . + b else . end
end
end;
 
def chinese_remainder(mods; remainders):
(reduce mods[] as $i (1; . * $i)) as $prod
| reduce range(0; mods|length) as $i
(0;
($prod/mods[$i]) as $p
| mul_inv($p; mods[$i]) as $mi
| if $mi == null then error("nogo: p=\($p) mods[\($i)]=\(mods[$i])")
else . + (remainders[$i] * $mi * $p)
end )
| . % $prod ;

Examples:

chinese_remainder([3,5,7]; [2,3,2])   
# => 23
chinese_remainder([100,23]; [19,0])
# => 1219
chinese_remainder([10,4,9]; [11,22,19])
# jq: error: nogo: p=36 mods[0]=10

Maple[edit]

This is a Maple built-in procedure, so it is trivial:

> chrem( [2, 3, 2], [3, 5, 7] );
23
 


Mathematica / Wolfram Language[edit]

Very easy, because it is a built-in function:

ChineseRemainder[{2, 3, 2}, {3, 5, 7}]
23

Nim[edit]

Translation of: C
proc mulInv(a0, b0): int =
var (a, b, x0) = (a0, b0, 0)
result = 1
if b == 1: return
while a > 1:
let q = a div b
a = a mod b
swap a, b
result = result - q * x0
swap x0, result
if result < 0: result += b0
 
proc chineseRemainder[T](n, a: T): int =
var prod = 1
var sum = 0
for x in n: prod *= x
 
for i in 0 .. <n.len:
let p = prod div n[i]
sum += a[i] * mulInv(p, n[i]) * p
 
sum mod prod
 
echo chineseRemainder([3,5,7], [2,3,2])

Output:

23

OCaml[edit]

This is using the Jane Street Ocaml Core library.

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')
 
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  

PARI/GP[edit]

chivec(residues, moduli)={
my(m=Mod(0,1));
for(i=1,#residues,
m=chinese(Mod(residues[i],moduli[i]),m)
);
lift(m)
};
chivec([2,3,2], [3,5,7])
Output:
23

Pari's chinese function takes a vector in the form [Mod(a1,n1), Mod(a2, n2), ...], so we can do this directly:

lift( chinese([Mod(2,3),Mod(3,5),Mod(2,7)]) )

or to take the residue/moduli array as above:

chivec(residues,moduli)={
lift(chinese(vector(#residues,i,Mod(residues[i],moduli[i]))))
}

Perl[edit]

There are at least three CPAN modules for this: ntheory (Math::Prime::Util), Math::ModInt, and Math::Pari. All three handle bigints.

Library: ntheory
use ntheory qw/chinese/;
say chinese([2,3], [3,5], [2,7]);
Output:
23

The function returns undef if no common residue class exists. The combined modulus can be obtained using the lcm function applied to the moduli (e.g. lcm(3,5,7) = 105 in the example above).

use Math::ModInt qw(mod);
use Math::ModInt::ChineseRemainder qw(cr_combine);
say cr_combine(mod(2,3),mod(3,5),mod(2,7));
Output:
mod(23, 105)

This returns a Math::ModInt object, which if no common residue class exists will be a special undefined object. The modulus and residue methods may be used to extract the integer components.

Non-pairwise-coprime[edit]

All three modules will also handle cases where the moduli are not pairwise co-prime but a solution exists, e.g.:

use ntheory qw/chinese/;
say chinese( [2328,16256], [410,5418] ), " mod ", lcm(16256,5418);
Output:
28450328 mod 44037504

Perl 6[edit]

Translation of: C
Works with: Rakudo version 2015.12
# 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);
Output:
23

PicoLisp[edit]

(de modinv (A B)
(let (B0 B X0 0 X1 1 Q 0 T1 0)
(while (< 1 A)
(setq
Q (/ A B)
T1 B
B (% A B)
A T1
T1 X0
X0 (- X1 (* Q X0))
X1 T1 ) )
(if (lt0 X1) (+ X1 B0) X1) ) )
 
(de chinrem (N A)
(let P (apply * N)
(%
(sum
'((N A)
(setq T1 (/ P N))
(* A (modinv T1 N) T1) )
N
A )
P ) ) )
 
(println
(chinrem (3 5 7) (2 3 2))
(chinrem (11 12 13) (10 4 12)) )
 
(bye)

Python[edit]

# Python 2.7
def chinese_remainder(n, a):
sum = 0
prod = reduce(lambda a, b: a*b, n)
 
for n_i, a_i in zip(n, a):
p = prod / n_i
sum += a_i * mul_inv(p, n_i) * p
return sum % prod
 
 
def mul_inv(a, b):
b0 = b
x0, x1 = 0, 1
if b == 1: return 1
while a > 1:
q = a / b
a, b = b, a%b
x0, x1 = x1 - q * x0, x0
if x1 < 0: x1 += b0
return x1
 
if __name__ == '__main__':
n = [3, 5, 7]
a = [2, 3, 2]
print chinese_remainder(n, a)
Output:
23

R[edit]

Translation of: C
mul_inv <- function(a, b)
{
b0 <- b
x0 <- 0L
x1 <- 1L
 
if (b == 1) return(1L)
while(a > 1){
q <- as.integer(a/b)
 
t <- b
b <- a %% b
a <- t
 
t <- x0
x0 <- x1 - q*x0
x1 <- t
}
 
if (x1 < 0) x1 <- x1 + b0
return(x1)
}
 
chinese_remainder <- function(n, a)
{
len <- length(n)
 
prod <- 1L
sum <- 0L
 
for (i in 1:len) prod <- prod * n[i]
 
for (i in 1:len){
p <- as.integer(prod / n[i])
sum <- sum + a[i] * mul_inv(p, n[i]) * p
}
 
return(sum %% prod)
}
 
n <- c(3L, 5L, 7L)
a <- c(2L, 3L, 2L)
 
chinese_remainder(n, a)
Output:
23


Racket[edit]

This is more of a demonstration of the built-in function "solve-chinese", than anything. A bit cheeky, I know... but if you've got a dog, why bark yourself?

Take a look in the "math/number-theory" package it's full of goodies! URL removed -- I can't be doing the Dutch recaptchas I'm getting.

#lang racket
(require (only-in math/number-theory solve-chinese))
(define as '(2 3 2))
(define ns '(3 5 7))
(solve-chinese as ns)
Output:
23

REXX[edit]

algebraic[edit]

/*REXX program demonstrates  Sun Tzu's  (or Sunzi's)  Chinese Remainder  Theorem.       */
parse arg Ns As . /*get optional arguments from the C.L. */
if Ns=='' | Ns=="," then Ns = '3,5,7' /*Ns not specified? Then use default.*/
if As=='' | As=="," then As = '2,3,2' /*As " " " " " */
say 'Ns: ' Ns
say 'As: ' As; say
Ns=space(translate(Ns, , ',')); #=words(Ns) /*elide any superfluous blanks from N's*/
As=space(translate(As, , ',')); _=words(As) /* " " " " " A's*/
if #\==_ then do; say "size of number sets don't match."; exit 131; end
if #==0 then do; say "size of the N set isn't valid."; exit 132; end
if _==0 then do; say "size of the A set isn't valid."; exit 133; end
N=1 /*the product─to─be for prod(n.j). */
do j=1 for # /*process each number for As and Ns. */
n.j=word(Ns,j); N=N*n.j /*get an N.j and calculate product. */
a.j=word(As,j) /* " " A.j from the As list. */
end /*j*/
 
do x=1 for N /*use a simple algebraic method. */
do i=1 for # /*process each N.i and A.i number.*/
if x//n.i\==a.i then iterate x /*is modulus correct for the number X ?*/
end /*i*/ /* [↑] limit solution to the product. */
say 'found a solution with X=' x /*display one possible solution. */
exit /*stick a fork in it, we're all done. */
end /*x*/
 
say 'no solution found.' /*oops, announce that solution ¬ found.*/
 

output

Ns:  3,5,7
As:  2,3,2

found a solution with X= 23

congruences sets[edit]

/*REXX program demonstrates  Sun Tzu's  (or Sunzi's)  Chinese Remainder  Theorem.       */
parse arg Ns As . /*get optional arguments from the C.L. */
if Ns=='' | Ns=="," then Ns = '3,5,7' /*Ns not specified? Then use default.*/
if As=='' | As=="," then As = '2,3,2' /*As " " " " " */
say 'Ns: ' Ns
say 'As: ' As; say
Ns=space(translate(Ns, , ',')); #=words(Ns) /*elide any superfluous blanks from N's*/
As=space(translate(As, , ',')); _=words(As) /* " " " " " A's*/
if #\==_ then do; say "size of number sets don't match."; exit 131; end
if #==0 then do; say "size of the N set isn't valid."; exit 132; end
if _==0 then do; say "size of the A set isn't valid."; exit 133; end
N=1 /*the product─to─be for prod(n.j). */
do j=1 for # /*process each number for As and Ns. */
n.j=word(Ns,j); N=N*n.j /*get an N.j and calculate product. */
a.j=word(As,j) /* " " A.j from the As list. */
end /*j*/
@.= /* [↓] converts congruences ───► sets.*/
do i=1 for #; _=a.i; @.i._=a.i; p=a.i
do N; p=p+n.i; @.i.p=p; end /*build a (array) list of modulo values*/
end /*i*/
/* [↓] find common number in the sets.*/
do x=1 for N; if @.1.x=='' then iterate /*locate a number. */
do v=2 to #; if @.v.x=='' then iterate x; end /*Is in all sets ? */
say 'found a solution with X=' x /*display one possible solution. */
exit /*stick a fork in it, we're all done. */
end /*x*/
 
say 'no solution found.' /*oops, announce that solution ¬ found.*/

output   is the same as the 1st REXX version.

Ruby[edit]

Brute-force.

 
def chinese_remainder(mods, remainders)
max = mods.inject( :* )
series = remainders.zip( mods ).map{|r,m| r.step( max, m ).to_a }
series.inject( :& ).first #returns nil when empty
end
 
p chinese_remainder([3,5,7], [2,3,2]) #=> 23
p chinese_remainder([10,4,9], [11,22,19]) #=> nil
 

Similar to above, but working with large(r) numbers.

 
def extended_gcd(a, b)
last_remainder, remainder = a.abs, b.abs
x, last_x, y, last_y = 0, 1, 1, 0
while remainder != 0
last_remainder, (quotient, remainder) = remainder, last_remainder.divmod(remainder)
x, last_x = last_x - quotient*x, x
y, last_y = last_y - quotient*y, y
end
return last_remainder, last_x * (a < 0 ? -1 : 1)
end
 
def invmod(e, et)
g, x = extended_gcd(e, et)
if g != 1
raise 'Multiplicative inverse modulo does not exist!'
end
x % et
end
 
def chinese_remainder(mods, remainders)
max = mods.inject( :* ) # product of all moduli
series = remainders.zip(mods).map{ |r,m| (r * max * invmod(max/m, m) / m) }
series.inject( :+ ) % max
end
 
p chinese_remainder([3,5,7], [2,3,2]) #=> 23
p chinese_remainder([17353461355013928499, 3882485124428619605195281, 13563122655762143587], [7631415079307304117, 1248561880341424820456626, 2756437267211517231]) #=> 937307771161836294247413550632295202816
p chinese_remainder([10,4,9], [11,22,19]) #=> nil
 

Seed7[edit]

$ include "seed7_05.s7i";
include "bigint.s7i";
 
const func integer: modInverse (in integer: a, in integer: b) is
return ord(modInverse(bigInteger conv a, bigInteger conv b));
 
const proc: main is func
local
const array integer: n is [] (3, 5, 7);
const array integer: a is [] (2, 3, 2);
var integer: num is 0;
var integer: prod is 1;
var integer: sum is 0;
var integer: index is 0;
begin
for num range n do
prod *:= num;
end for;
for key index range a do
num := prod div n[index];
sum +:= a[index] * modInverse(num, n[index]) * num;
end for;
writeln(sum mod prod);
end func;
Output:
23

Sidef[edit]

Translation of: Perl 6
func mul_inv(a, b) {
b == 1 && return 1;
var (b0, x0, x1) = (0, 0, 1);
while (a > 1) {
(a, b, x0, x1) = (b, a % b, x1 - x0*int(a / b), x0);
};
x1 < 0 ? x1+b0 : x1;
}
 
func chinese_remainder(*n) {
var N = n«*»;
func (*a) {
n.range.map { |i|
var p = int(N / n[i]);
a[i] * mul_inv(p, n[i]) * p;
}.sum
}
}
 
say chinese_remainder(3, 5, 7)(2, 3, 2);
Output:
23

Tcl[edit]

Translation of: C
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}]
Output:
23

uBasic/4tH[edit]

Translation of: C
@(000) = 3 : @(001) = 5 : @(002) = 7
@(100) = 2 : @(101) = 3 : @(102) = 2
 
Print Func (_Chinese_Remainder (3))
 
' -------------------------------------
 
@(000) = 11 : @(001) = 12 : @(002) = 13
@(100) = 10 : @(101) = 04 : @(102) = 12
 
Print Func (_Chinese_Remainder (3))
 
' -------------------------------------
 
End
 
' returns x where (a * x) % b == 1
_Mul_Inv Param (2) ' ( a b -- n)
Local (4)
 
c@ = b@
d@ = 0
e@ = 1
 
If b@ = 1 Then Return (1)
 
Do While a@ > 1
f@ = a@ / b@
Push b@ : b@ = a@ % b@ : a@ = Pop()
Push d@ : d@ = e@ - f@ * d@ : e@ = Pop()
Loop
 
If e@ < 0 Then e@ = e@ + c@
 
Return (e@)
 
 
_Chinese_Remainder Param (1) ' ( len -- n)
Local (5)
 
b@ = 1
c@ = 0
 
For d@ = 0 Step 1 While d@ < a@
b@ = b@ * @(d@)
Next
 
For d@ = 0 Step 1 While d@ < a@
e@ = b@ / @(d@)
c@ = c@ + (@(100 + d@) * Func (_Mul_Inv (e@, @(d@))) * e@)
Next
 
Return (c@ % b@)
 
Output:
23
1000

0 OK, 0:1034

zkl[edit]

Translation of: Go

Using the GMP library, gcdExt is the Extended Euclidean algorithm.

var BN=Import("zklBigNum"), one=BN(1);
 
fcn crt(xs,ys){
p:=xs.reduce('*,BN(1));
X:=BN(0);
foreach x,y in (xs.zip(ys)){
q:=p/x;
z,s,_:=q.gcdExt(x);
if(z!=one) throw(Exception.ValueError("%d not coprime".fmt(x)));
X+=y*s*q;
}
return(X % p);
}
println(crt(T(3,5,7),  T(2,3,2)));    //-->23
println(crt(T(11,12,13),T(10,4,12))); //-->1000
println(crt(T(11,22,19), T(10,4,9))); //-->ValueError: 11 not coprime

ZX Spectrum Basic[edit]

Translation of: C
10 DIM n(3): DIM a(3)
20 FOR i=1 TO 3
30 READ n(i),a(i)
40 NEXT i
50 DATA 3,2,5,3,7,2
100 LET prod=1: LET sum=0
110 FOR i=1 TO 3: LET prod=prod*n(i): NEXT i
120 FOR i=1 TO 3
130 LET p=INT (prod/n(i)): LET a=p: LET b=n(i)
140 GO SUB 1000
150 LET sum=sum+a(i)*x1*p
160 NEXT i
170 PRINT FN m(sum,prod)
180 STOP
200 DEF FN m(a,b)=a-INT (a/b)*b: REM Modulus function
1000 LET b0=b: LET x0=0: LET x1=1
1010 IF b=1 THEN RETURN
1020 IF a<=1 THEN GO TO 1100
1030 LET q=INT (a/b)
1040 LET t=b: LET b=FN m(a,b): LET a=t
1050 LET t=x0: LET x0=x1-q*x0: LET x1=t
1060 GO TO 1020
1100 IF x1<0 THEN LET x1=x1+b0
1110 RETURN
23