Chinese remainder theorem: Difference between revisions

Content deleted Content added

Inline

Revision as of 21:35, 23 December 2014

Suppose Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_1} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_2} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \ldots} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_k} are positive integers that are pairwise coprime. Then, for any given sequence of integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_1} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_2} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \dots} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_k} , there exists an integer Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x} solving the following system of simultaneous congruences.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \begin{align} x &\equiv a_1 \pmod{n_1} \\ x &\equiv a_2 \pmod{n_2} \\ &{}\ \ \vdots \\ x &\equiv a_k \pmod{n_k} \end{align}}

Furthermore, all solutions Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x} of this system are congruent modulo the product, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N=n_1n_2\ldots 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s} where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0 \leq s \leq n_1n_2\ldots n_k} .

Show the functionality of this program by printing the result such that the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n} 's are Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [3,5,7]} and the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} 's are Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [2,3,2]} .

Algorithm: The following algorithm only applies if the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_i} 's are pairwise coprime.

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x \equiv a_i \pmod{n_i} \quad\mathrm{for}\; i = 1, \ldots, k}

Again, to begin, the product Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N = n_1n_2 \ldots n_k} is defined. Then a solution Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x} can be found as follows.

For each Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle i} , the integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n_i} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N/n_i} are coprime. Using the Extended Euclidean algorithm we can find integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_i} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle s_i} such that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle r_i n_i + s_i N/n_i = 1} . Then, one solution to the system of simultaneous congruences is:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x = \sum_{i=1}^k a_i s_i N/n_i}

and the minimal solution,

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x \pmod{N}} .

Ada

Using the package Mod_Inv from [[1]].

<lang Ada>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;</lang>

Bracmat

Translation of: C

<lang bracmat>( ( 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)) );</lang> Output:

C

When n are not pairwise 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>

Common Lisp

Using function invmod from [[2]]. <lang lisp> (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)))))

</lang>

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

Translation of: Python

<lang d>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;

}</lang>

Output:

Erlang

Translation of: OCaml

<lang erlang>-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.</lang>

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

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

Go

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. <lang go>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))

}</lang>

Output:

Two values, the solution x and an error value.

23 <nil>

Haskell

Translation of: Erlang

<lang haskell>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

</lang>

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

Translation of: Python

with error check added.

Works in both languages: <lang unicon>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</lang>

Output:

->crt
23
No solution!
->

J

Solution (brute force):<lang j> crt =: (1 + ] - {:@:[ -: {.@:[ | ])^:_&0@:,:</lang> Example: <lang j> 3 5 7 crt 2 3 2 23

  11 12 13 crt 10 4 12

1000</lang> 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.

jq

This implementation is similar to the one in C, but raises an error if there is no solution, as illustrated in the last example. <lang jq># 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 ;</lang>

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

This is a Maple built-in procedure, so it is trivial: <lang Maple>> chrem( [2, 3, 2], [3, 5, 7] );

</lang>

Mathematica

Very easy, because it is a built-in function: <lang Mathematica >ChineseRemainder[{2, 3, 2}, {3, 5, 7}] 23</lang>

Nimrod

Translation of: C

<lang nimrod>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])</lang> Output:

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

PARI/GP

<lang parigp>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])</lang>

Output:

Pari's chinese function takes a vector in the form [Mod(a1,n1), Mod(a2, n2), ...], so we can do this directly: <lang parigp>lift( chinese([Mod(2,3),Mod(3,5),Mod(2,7)]) )</lang> or to take the residue/moduli array as above: <lang parigp>chivec(residues,moduli)={

 lift(chinese(vector(#residues,i,Mod(residues[i],moduli[i]))))

}</lang>

Perl

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

Library: ntheory

<lang perl>use ntheory qw/chinese/; say chinese([2,3], [3,5], [2,7]);</lang>

Output:

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

<lang perl>use Math::ModInt qw(mod); use Math::ModInt::ChineseRemainder qw(cr_combine); say cr_combine(mod(2,3),mod(3,5),mod(2,7));</lang>

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

All three modules will also handle cases where the moduli are not pairwise co-prime but a solution exists, e.g.: <lang perl>use ntheory qw/chinese/; say chinese( [2328,16256], [410,5418] ), " mod ", lcm(16256,5418);</lang>

Output:

28450328 mod 44037504

Perl 6

Translation of: C

<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:

PicoLisp

<lang PicoLisp>(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)</lang>

Python

Translation of: C

<lang python>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

def chinese_remainder(n, a, lena):

   p = i = prod = 1; sm = 0
   for i in range(lena): prod *= n[i]
   for i in range(lena):
       p = prod / n[i]
       sm += a[i] * mul_inv(p, n[i]) * p
   return sm % prod

if __name__ == '__main__':

   n = [3, 5, 7]
   a = [2, 3, 2]
   print(chinese_remainder(n, a, len(n)))</lang>

Output:

R

Translation of: C

<lang rsplus>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)</lang>

Output:

Racket

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>#lang racket (require (only-in math/number-theory solve-chinese)) (define as '(2 3 2)) (define ns '(3 5 7)) (solve-chinese as ns)</lang>

Output:

REXX

algebraic

<lang rexx>/*REXX program uses the Chinese Remainder Theorem (Sun Tzu). */ parse arg Ns As . /*get optional arguments from CL.*/ if Ns== then Ns = '3,5,7' /*Ns not specified? Use default.*/ if As== then As = '2,3,2' /*As " " " " */ Ns=space(translate(Ns,,',')); #=words(Ns) /*elide superfluous blanks*/ As=space(translate(As,,',')); _=words(As) /* " " " */ 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, Ns.*/
     n.j=word(Ns,j);  N=N*n.j         /*get an  N.j  and calculate prod*/
     a.j=word(As,j)                   /* "   "  A.j  from the  As.     */
     end   /*j*/

 do    x=1  for N                     /*use a simple algebraic method. */
    do i=1  for #                     /*process each   A.i   number.   */
    if x//n.i\==a.i  then iterate x   /*is the modulus correct for F ? */
    end   /*i*/                       /* [↑]  limit solution to product*/
 say 'found a solution with x=' x     /*announce a possible solution.  */
 exit                                 /*stick a fork in it, we're done.*/
 end   /*x*/
                                      /*stick a fork in it, we're done.*/

say 'no solution found.' /*oops, announce that ¬ found. */</lang> output

found a solution with x= 23

congruences sets

<lang rexx>/*REXX program uses the Chinese Remainder Theorem (Sun Tzu). */ parse arg Ns As . /*get optional arguments from CL.*/ if Ns== then Ns = '3,5,7' /*Ns not specified? Use default.*/ if As== then As = '2,3,2' /*As " " " " */ Ns=space(translate(Ns,,',')); #=words(Ns) /*elide superfluous blanks*/ As=space(translate(As,,',')); _=words(As) /* " " " */ 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, Ns.*/
     n.j=word(Ns,j);  N=N*n.j         /*get an  N.j  and calculate prod*/
     a.j=word(As,j)                   /* "   "  A.j  from the  As.     */
     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 list of modulo values. */
   end   /*i*/
                                      /* [↓] find common number in sets*/
 do   x=1  for N;  if @.1.x==  then iterate          /*find a number.*/
   do v=2  to #;   if @.v.x==  then iterate x;  end  /*In all sets ? */
 say 'found a solution with X=' x     /*we found the lowest solution.  */
 exit                                 /*stick a fork in it, we're done.*/
 end   /*x*/
                                      /*stick a fork in it, we're done.*/

say 'no solution found.' /*oops, there's not a solution. */</lang> output

found a solution with X= 23

Ruby

Brute-force. <lang ruby> 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 </lang>

Seed7

<lang seed7>$ 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;</lang>

Output:

Tcl

Translation of: C

<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:

uBasic/4tH

Translation of: C

<lang>@(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@) </lang>

Output:

23
1000

0 OK, 0:1034

zkl

Translation of: Go

Using the GMP library, gcdExt is the Extended Euclidean algorithm. <lang zkl>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);

}</lang> <lang zkl>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</lang>