Chinese remainder theorem: Difference between revisions

Content added Content deleted

Inline

Revision as of 07:06, 9 August 2023

Suppose $n_{1}$ , $n_{2}$ , $\ldots$ , $n_{k}$ are positive integers that are pairwise co-prime.

Then, for any given sequence of integers $a_{1}$ , $a_{2}$ , $\dots$ , $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_{1}n_{2}\ldots n_{k}$ .

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 $s$ where $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$ 's are $[3,5,7]$ and the $a$ 's are $[2,3,2]$ .

Algorithm: The following algorithm only applies if the $n_{i}$ 's are pairwise co-prime.

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 $N=n_{1}n_{2}\ldots n_{k}$ is defined.

Then a solution $x$ can be found as follows:

For each $i$ , the integers $n_{i}$ and $N/n_{i}$ are co-prime.

Using the Extended Euclidean algorithm, we can find integers $r_{i}$ and $s_{i}$ such that $r_{i}n_{i}+s_{i}N/n_{i}=1$ .

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

x=\sum _{i=1}^{k}a_{i}s_{i}N/n_{i}

and the minimal solution,

x{\pmod {N}}

.

360 Assembly

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

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits
or android 64 bits with application Termux

/* ARM assembly AARCH64 Raspberry PI 3B */
/*  program chineserem64.s   */

/************************************/
/* Constantes                       */
/************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc" 


/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessResult:       .asciz "Result = "
szCarriageReturn:   .asciz "\n"
.align 2
arrayN:           .quad 3,5,7
arrayA:           .quad 2,3,2
      .equ ARRAYSIZE,  (. - arrayA)/8
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss  
sZoneConv:           .skip 24
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:
    ldr x0,qAdrarrayN            // N array address
    ldr x1,qAdrarrayA            // A array address
    mov x2,#ARRAYSIZE            //  array size
    bl chineseremainder
 
    ldr x1,qAdrsZoneConv       
    bl conversion10             // call décimal conversion
    mov x0,#3
    ldr x1,qAdrszMessResult
    ldr x2,qAdrsZoneConv        // insert conversion in message
    ldr x3,qAdrszCarriageReturn
    bl displayStrings           // display message

100:                            // standard end of the program 
    mov x0, #0                  // return code
    mov x8,EXIT 
    svc #0                      // perform the system call
qAdrszCarriageReturn:    .quad szCarriageReturn
qAdrsZoneConv:           .quad sZoneConv  
qAdrszMessResult:        .quad szMessResult
qAdrarrayA:              .quad arrayA
qAdrarrayN:              .quad arrayN

/******************************************************************/
/*     compute chinese remainder                                  */ 
/******************************************************************/
/* x0 contains n array address */
/* x1 contains a array address */
/* x2 contains array size      */
chineseremainder:
    stp x1,lr,[sp,-16]!       // save  registers 
    stp x2,x3,[sp,-16]!       // save  registers 
    stp x4,x5,[sp,-16]!       // save  registers 
    stp x6,x7,[sp,-16]!       // save  registers 
    stp x8,x9,[sp,-16]!       // save  registers 
    mov x4,#1                 // product
    mov x5,#0                 // sum
    mov x6,#0                 // indice
1: 
    ldr x3,[x0,x6,lsl #3]     // load a value 
    mul x4,x3,x4              // compute product
    add x6,x6,#1
    cmp x6,x2
    blt 1b
    mov x6,#0
    mov x7,x0                 // save entry
    mov x8,x1
    mov x9,x2
2: 
    mov x0,x4                 // product
    ldr x1,[x7,x6,lsl #3]     // value of n
    sdiv x2,x0,x1
    mov x0,x2                 //  p
    bl inverseModulo
    mul x0,x2,x0              // = product / n * invmod
    ldr x3,[x8,x6,lsl #3]     //  value a
    madd x5,x0,x3,x5          // sum = sum + (result1 * a)
    add x6,x6,#1
    cmp x6,x9
    blt 2b
    sdiv x1,x5,x4             // divide sum by produc
    msub x0,x1,x4,x5          // compute remainder
100:
    ldp x8,x9,[sp],16        // restaur  registers 
    ldp x6,x7,[sp],16        // restaur  registers 
    ldp x4,x5,[sp],16        // restaur  registers 
    ldp x2,x3,[sp],16        // restaur  registers 
    ldp x1,lr,[sp],16            // restaur  registers
    ret 
/***************************************************/
/*   Calcul modulo inverse                            */
/***************************************************/
/* x0 cont.quad number, x1 modulo              */
/* x0 return result               */
inverseModulo:
    stp x1,lr,[sp,-16]!          // save  registers 
    stp x2,x3,[sp,-16]!          // save  registers 
    stp x4,x5,[sp,-16]!          // save  registers 
    stp x6,x7,[sp,-16]!          // save  registers 
    mov x7,x1            // save Modulo
    mov x6,x1            // A   x0=B
    mov x4,#1            // X
    mov x5,#0            // Y
1:   // 
    cmp x0,#0            // B = 0
    beq 2f
    mov x1,x0            // T = B
    mov x0,x6            // A
    sdiv x2,x0,x1        // A / T
    msub x0,x2,x1,x0     // B and x2=Q
    mov x6,x1            // A=T
    mov x1,x4            // T=X
    msub x4,x2,x1,x5     // X=Y-(Q*T)
    mov x5,x1            // Y=T
    b 1b
2:
    add x7,x7,x5         // = Y + N
    cmp x5,#0            // Y > 0
    bge 3f
    mov x0,x7
    b 100f
3:
    mov x0,x5
100:
    ldp x6,x7,[sp],16        // restaur  registers 
    ldp x4,x5,[sp],16        // restaur  registers 
    ldp x2,x3,[sp],16        // restaur  registers 
    ldp x1,lr,[sp],16            // restaur  registers
    ret 
/***************************************************/
/*   display multi strings                    */
/***************************************************/
/* x0  contains number strings address */
/* x1 address string1 */
/* x2 address string2 */
/* x3 address string3 */
/* other address on the stack */
/* thinck to add  number other address * 4 to add to the stack */
displayStrings:            // INFO:  affichageStrings
    stp x1,lr,[sp,-16]!    // save  registers 
    stp x2,x3,[sp,-16]!    // save  registers 
    stp x4,x5,[sp,-16]!    // save  registers 
    add fp,sp,#48          // save paraméters address (6 registers saved * 8 bytes)
    mov x4,x0              // save strings number
    cmp x4,#0              // 0 string -> end
    ble 100f
    mov x0,x1              // string 1
    bl affichageMess
    cmp x4,#1              // number > 1
    ble 100f
    mov x0,x2
    bl affichageMess
    cmp x4,#2
    ble 100f
    mov x0,x3
    bl affichageMess
    cmp x4,#3
    ble 100f
    mov x3,#3
    sub x2,x4,#4
1:                          // loop extract address string on stack
    ldr x0,[fp,x2,lsl #3]
    bl affichageMess
    subs x2,x2,#1
    bge 1b
100:
    ldp x4,x5,[sp],16        // restaur  registers 
    ldp x2,x3,[sp],16        // restaur  registers 
    ldp x1,lr,[sp],16        // restaur  registers
    ret 
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
/* for this file see task include a file in language AArch64 assembly */
.include "../includeARM64.inc"

Output:

Result = 23

Action!

INT FUNC MulInv(INT a,b)
  INT b0,x0,x1,q,tmp

  IF b=1 THEN RETURN (1) FI

  b0=b x0=0 x1=1
  WHILE a>1
  DO
    q=a/b

    tmp=b
    b=a MOD b
    a=tmp
    
    tmp=x0
    x0=x1-q*x0
    x1=tmp
  OD
  IF x1<0 THEN
    x1==+b0
  FI
RETURN (x1)

INT FUNC ChineseRemainder(BYTE ARRAY n,a BYTE len)
  INT prod,sum,p,m
  BYTE i

  prod=1 sum=0
  FOR i=0 TO len-1
  DO
    prod==*n(i)
  OD
  FOR i=0 TO len-1
  DO
    p=prod/n(i)
    m=MulInv(p,n(i))
    sum==+a(i)*m*p
  OD
RETURN (sum MOD prod)

PROC Main()
  BYTE ARRAY n=[3 5 7],a=[2 3 2]
  INT res

  res=ChineseRemainder(n,a,3)
  PrintI(res)
RETURN

Output:

Screenshot from Atari 8-bit computer

Ada

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;

ARM Assembly

Works with: as version Raspberry Pi
or android 32 bits with application Termux

/* ARM assembly Raspberry PI or android with termux */
/*  program chineserem.s   */

 /* REMARK 1 : this program use routines in a include file 
   see task Include a file language arm assembly 
   for the routine affichageMess conversion10 
   see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes                       */
/************************************/
.include "../constantes.inc"


/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessResult:       .asciz "Result = "
szCarriageReturn:   .asciz "\n"
.align 2
arrayN:           .int 3,5,7
arrayA:           .int 2,3,2
      .equ ARRAYSIZE,  (. - arrayA)/4
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss  
sZoneConv:           .skip 24
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:
    ldr r0,iAdrarrayN            @ N array address
    ldr r1,iAdrarrayA            @ A array address
    mov r2,#ARRAYSIZE            @  array size
    bl chineseremainder
 
    ldr r1,iAdrsZoneConv       
    bl conversion10             @ call décimal conversion
    mov r0,#3
    ldr r1,iAdrszMessResult
    ldr r2,iAdrsZoneConv        @ insert conversion in message
    ldr r3,iAdrszCarriageReturn
    bl displayStrings           @ display message

100:                            @ standard end of the program 
    mov r0, #0                  @ return code
    mov r7, #EXIT               @ request to exit program
    svc #0                      @ perform the system call
iAdrszCarriageReturn:    .int szCarriageReturn
iAdrsZoneConv:           .int sZoneConv  
iAdrszMessResult:        .int szMessResult
iAdrarrayA:              .int arrayA
iAdrarrayN:              .int arrayN

/******************************************************************/
/*     compute chinese remainder                                  */ 
/******************************************************************/
/* r0 contains n array address */
/* r1 contains a array address */
/* r2 contains array size      */
chineseremainder:
    push {r1-r9,lr}           @ save  registers 
    mov r4,#1                 @ product
    mov r5,#0                 @ sum
    mov r6,#0                 @ indice
1: 
    ldr r3,[r0,r6,lsl #2]     @ load a value 
    mul r4,r3,r4              @ compute product
    add r6,#1
    cmp r6,r2
    blt 1b
    mov r6,#0
    mov r7,r0                 @ save entry
    mov r8,r1
    mov r9,r2
2: 
    mov r0,r4                 @ product
    ldr r1,[r7,r6,lsl #2]     @ value of n
    bl division
    mov r0,r2                 @  p
    bl inverseModulo
    mul r0,r2,r0              @ = product / n * invmod
    ldr r3,[r8,r6,lsl #2]     @  value a
    mla r5,r0,r3,r5           @ sum = sum + (result1 * a)
    add r6,#1
    cmp r6,r9
    blt 2b
    mov r0,r5                  @ sum
    mov r1,r4                  @ product
    bl division
    mov r0,r3
    
100:
    pop {r1-r9,pc}              @ restaur registers
/***************************************************/
/*   Calcul modulo inverse                            */
/***************************************************/
/* r0 containt number, r1 modulo              */
/* x0 return result               */
inverseModulo:
    push {r1-r7,lr}       @ save  registers 
    mov r7,r1            // save Modulo
    mov r6,r1            // A   r0=B
    mov r4,#1            // X
    mov r5,#0            // Y
1:   // 
    cmp r0,#0            // B = 0
    beq 2f
    mov r1,r0            // T = B
    mov r0,r6            // A
    bl division          // A / T
    mov r0,r3            // B and r2=Q
    mov r6,r1            // A=T
    mov r1,r4            // T=X
    mls r4,r2,r1,r5      // X=Y-(Q*T)
    mov r5,r1            // Y=T
    b 1b
2:
    add r7,r7,r5         // = Y + N
    cmp r5,#0            // Y > 0
    bge 3f
    mov r0,r7
    b 100f
3:
    mov r0,r5
100:
   pop {r1-r7,pc}
/***************************************************/
/*   display multi strings                    */
/***************************************************/
/* r0  contains number strings address */
/* r1 address string1 */
/* r2 address string2 */
/* r3 address string3 */
/* other address on the stack */
/* thinck to add  number other address * 4 to add to the stack */
displayStrings:            @ INFO:  affichageStrings
    push {r1-r4,fp,lr}     @ save des registres
    add fp,sp,#24          @ save paraméters address (6 registers saved * 4 bytes)
    mov r4,r0              @ save strings number
    cmp r4,#0              @ 0 string -> end
    ble 100f
    mov r0,r1              @ string 1
    bl affichageMess
    cmp r4,#1              @ number > 1
    ble 100f
    mov r0,r2
    bl affichageMess
    cmp r4,#2
    ble 100f
    mov r0,r3
    bl affichageMess
    cmp r4,#3
    ble 100f
    mov r3,#3
    sub r2,r4,#4
1:                            @ loop extract address string on stack
    ldr r0,[fp,r2,lsl #2]
    bl affichageMess
    subs r2,#1
    bge 1b
100:
    pop {r1-r4,fp,pc}
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
.include "../affichage.inc"

Output:

Result = 23

BQN

Multiplicative Modular inverse function taken from BQNcrate.

MulInv←⊣|·⊑{0=𝕨?1‿0;(⌽-(0⋈𝕩⌊∘÷𝕨)⊸×)𝕨𝕊˜𝕨|𝕩}
ChRem←{
  num 𝕊 rem:
  prod←×´num
  prod|+´rem×(⊢×num⊸(MulInv¨))prod⌊∘÷num
}

•Show 3‿5‿7 ChRem 2‿3‿2
•Show 10‿4‿9 ChRem 11‿22‿19

23
172

C

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

C#

using System;
using System.Linq;

namespace ChineseRemainderTheorem
{
    class Program
    {
        static void Main(string[] args)
        {
            int[] n = { 3, 5, 7 };
            int[] a = { 2, 3, 2 };

            int result = ChineseRemainderTheorem.Solve(n, a);

            int counter = 0;
            int maxCount = n.Length - 1;
            while (counter <= maxCount)
            {
                Console.WriteLine($"{result} ≡ {a[counter]} (mod {n[counter]})");
                counter++;
            }
        }
    }

    public static class ChineseRemainderTheorem
    {
        public static int Solve(int[] n, int[] a)
        {
            int prod = n.Aggregate(1, (i, j) => i * j);
            int p;
            int sm = 0;
            for (int i = 0; i < n.Length; i++)
            {
                p = prod / n[i];
                sm += a[i] * ModularMultiplicativeInverse(p, n[i]) * p;
            }
            return sm % prod;
        }

        private static int ModularMultiplicativeInverse(int a, int mod)
        {
            int b = a % mod;
            for (int x = 1; x < mod; x++)
            {
                if ((b * x) % mod == 1)
                {
                    return x;
                }
            }
            return 1;
        }
    }
}

Common Lisp

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]

Crystal

Translation of: Ruby

def extended_gcd(a, b)
    last_remainder, remainder = a.abs, b.abs
    x, last_x = 0, 1

    until remainder == 0
        tmp = remainder
        quotient, remainder = last_remainder.divmod(remainder)
        last_remainder = tmp
        x, last_x = last_x - quotient * x, x
    end

    return last_remainder, last_x * (a < 0 ? -1 : 1)
end


def invmod(e, et)
    g, x = extended_gcd(e, et)
    unless g == 1
        raise "Multiplicative inverse modulo does not exist"
    end
    return x % et
end


def chinese_remainder(mods, remainders)
    max = mods.product
    series = remainders.zip(mods).map { |r, m| r * max * invmod(max // m, m) // m }
    return series.sum % max
end


puts chinese_remainder([3, 5, 7], [2, 3, 2])
puts chinese_remainder([5, 7, 9, 11], [1, 2, 3, 4])

Output:

23
1731

EchoLisp

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

Translation of: Ruby

Works with: Elixir version 1.2

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) |> MapSet.new end)
    |> Enum.reduce(fn set,acc -> MapSet.intersection(set, acc) end)
    |> MapSet.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

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

F#

sieving

let rec sieve cs x N =
    match cs with
    | [] -> Some(x)
    | (a,n)::rest ->
        let arrProgress = Seq.unfold (fun x -> Some(x, x+N)) x
        let firstXmodNequalA = Seq.tryFind (fun x -> a = x % n)
        match firstXmodNequalA (Seq.take n arrProgress) with
        | None -> None
        | Some(x) -> sieve rest x (N*n)

[ [(2,3);(3,5);(2,7)];
  [(10,11); (4,22); (9,19)];
  [(10,11); (4,12); (12,13)] ]
|> List.iter (fun congruences ->
    let cs =
        congruences
        |> List.map (fun (a,n) -> (a % n, n))
        |> List.sortBy (snd>>(~-)) 
    let an = List.head cs
    match sieve (List.tail cs) (fst an) (snd an) with
    | None    -> printfn "no solution"
    | Some(x) -> printfn "result = %i" x
)

Output:

result = 23
no solution
result = 1000

Or for those who prefer unsieved

This uses Greatest_common_divisor#F.23 to verify valid input, can be simplified if you know input has a solution.
This uses Modular_inverse#F.23

//Chinese Division Theorem: Nigel Galloway: April 3rd., 2017
let CD n g =
  match Seq.fold(fun n g->if (gcd n g)=1 then n*g else 0) 1 g with
  |0 -> None
  |fN-> Some ((Seq.fold2(fun n i g -> n+i*(fN/g)*(MI g ((fN/g)%g))) 0 n g)%fN)

Output:

CD [10;4;12] [11;12;13] -> Some 1000
CD [10;4;9] [11;22;19]  -> None
CD [2;3;2] [3;5;7]      -> Some 23

Forth

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

Fortran

Written in Fortran-77 style (fixed form, GO TO), directly translated from the problem description.

* RC task: use the Chinese Remainder Theorem to solve a system of congruences.

      FUNCTION crt(n, residues, moduli)
      IMPLICIT INTEGER (A-Z)
      DIMENSION residues(n), moduli(n)

      p = product(moduli)
      crt = 0
      DO 10 i = 1, n
         m = p/moduli(i)
         CALL egcd(moduli(i), m, r, s, gcd)
         IF (gcd .ne. 1) GO TO 20 ! error exit
10       crt = crt + residues(i)*s*m
      crt = modulo(crt, p)
      RETURN

20    crt = -1 ! will never be negative, so flag an error
      END


* Compute egcd(a, b), returning x, y, g s.t.
*   g = gcd(a, b) and a*x + b*y = g
*   
      SUBROUTINE egcd(a, b, x, y, g)
      IMPLICIT INTEGER (A-Z)

      g = a
      u = 0
      v = 1
      w = b
      x = 1
      y = 0

1     IF (w .eq. 0) RETURN
      q = g/w
      u next = x - q*u
      v next = y - q*v
      w next = g - q*w
      x = u
      y = v
      g = w
      u = u next
      v = v next
      w = w next
      GO TO 1
      END


      PROGRAM Chinese Remainder
      IMPLICIT INTEGER (A-Z)

      PRINT *, crt(3, [2, 3, 2], [3, 5, 7])
      PRINT *, crt(3, [2, 3, 2], [3, 6, 7]) ! no solution

      END

Output:

          23
          -1

This example solves an extended version of the Chinese Remainder theorem by allowing an optional third parameter d which defaults to 0 and is an integer. The solution returned is the smallest solution >= d. (This optional parameter is common in many/most real-world applications of the Chinese Remainder Theorem.)

This program also works with arbitrarily-large integers and peforms efficiently due to Frink's built-in modInverse function.

Input is validated and useful error messages are emitted if the input data is invalid. If a solution cannot be found, this returns the special value undef.

/** arguments:
       [r, m, d=0]  where r and m are arrays of the remainder terms r and the
               modulus terms m respectively.  These must be of the same length.
    returns
       x, the unique solution mod N where N is the product of all the M terms where x &gt;= d.
*/
ChineseRemainder[r, m, d=0] :=
{
   if length[r] != length[m]
   {
      println["ChineseRemainder:  r and m must be arrays of the same length."]
      return undef
   }
   
   N = product[m]

   y = new array
   z = new array
   x = 0
   for i = rangeOf[m]
   {
      y@i = N / m@i
      z@i = modInverse[y@i, m@i]
      if z@i == undef
      {
         println["ChineseRemainder:  modInverse returned undef for modInverse[" + y@i + ", " + m@i + "]"]
         return undef
      }
      
      x = x + r@i y@i z@i
   }

   xp = x mod N
   f = d div N
   r = f * N + xp
   if r < d
      r = r + N

   return r
}

println[ChineseRemainder[[2,3,2],[3,5,7]] ]

Output:

FunL

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:

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.

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>

Groovy

Translation of: Java

class ChineseRemainderTheorem {
    static int chineseRemainder(int[] n, int[] a) {
        int prod = 1
        for (int i = 0; i < n.length; i++) {
            prod *= n[i]
        }

        int p, sm = 0
        for (int i = 0; i < n.length; i++) {
            p = prod.intdiv(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.intdiv(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
    }

    static void main(String[] args) {
        int[] n = [3, 5, 7]
        int[] a = [2, 3, 2]
        println(chineseRemainder(n, a))
    }
}

Output:

Haskell

Translation of: Erlang

import Control.Monad (zipWithM)

egcd :: Int -> Int -> (Int, Int)
egcd _ 0 = (1, 0)
egcd a b = (t, s - q * t)
  where
    (s, t) = egcd b r
    (q, r) = a `quotRem` b

modInv :: Int -> Int -> Either String Int
modInv a b =
  case egcd a b of
    (x, y)
      | a * x + b * y == 1 -> Right x
      | otherwise ->
        Left $ "No modular inverse for " ++ show a ++ " and " ++ show b

chineseRemainder :: [Int] -> [Int] -> Either String Int
chineseRemainder residues modulii =
  zipWithM modInv crtModulii modulii >>=
  (Right . (`mod` modPI) . sum . zipWith (*) crtModulii . zipWith (*) residues)
  where
    modPI = product modulii
    crtModulii = (modPI `div`) <$> modulii

main :: IO ()
main =
  mapM_ (putStrLn . either id show) $
  uncurry chineseRemainder <$>
  [ ([10, 4, 12], [11, 12, 13])
  , ([10, 4, 9], [11, 22, 19])
  , ([2, 3, 2], [3, 5, 7])
  ]

Output:

1000
No modular inverse for 418 and 11
23

Icon and Unicon

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

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

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

JavaScript

function crt(num, rem) {
  let sum = 0;
  const prod = num.reduce((a, c) => a * c, 1);

  for (let i = 0; i < num.length; i++) {
    const [ni, ri] = [num[i], rem[i]];
    const p = Math.floor(prod / ni);
    sum += ri * p * mulInv(p, ni);
  }
  return sum % prod;
}

function mulInv(a, b) {
  const b0 = b;
  let [x0, x1] = [0, 1];

  if (b === 1) {
    return 1;
  }
  while (a > 1) {
    const q = Math.floor(a / b);
    [a, b] = [b, a % b];
    [x0, x1] = [x1 - q * x0, x0];
  }
  if (x1 < 0) {
    x1 += b0;
  }
  return x1;
}

console.log(crt([3,5,7], [2,3,2]))

Output:

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.

# 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

Julia

Works with: Julia version 1.2

function chineseremainder(n::Array, a::Array)
    Π = prod(n)
    mod(sum(ai * invmod(Π ÷ ni, ni) * (Π ÷ ni) for (ni, ai) in zip(n, a)), Π)
end

@show chineseremainder([3, 5, 7], [2, 3, 2])

Output:

chineseremainder([3, 5, 7], [2, 3, 2]) = 23

Kotlin

Translation of: C

// version 1.1.2

/* returns x where (a * x) % b == 1 */
fun multInv(a: Int, b: Int): Int {
    if (b == 1) return 1
    var aa = a
    var bb = b
    var x0 = 0
    var x1 = 1
    while (aa > 1) {
        val q = aa / bb
        var t = bb
        bb = aa % bb
        aa = t
        t = x0
        x0 = x1 - q * x0
        x1 = t
    }
    if (x1 < 0) x1 += b
    return x1
} 

fun chineseRemainder(n: IntArray, a: IntArray): Int {
    val prod = n.fold(1) { acc, i -> acc * i }
    var sum = 0
    for (i in 0 until n.size) {
        val p = prod / n[i]
        sum += a[i] * multInv(p, n[i]) * p
    }
    return sum % prod
}

fun main(args: Array<String>) {
    val n = intArrayOf(3, 5, 7)
    val a = intArrayOf(2, 3, 2)
    println(chineseRemainder(n, a))
}

Output:

Lobster

import std

def extended_gcd(a, b):
    var s = 0
    var old_s = 1
    var t = 1
    var old_t = 0
    var r = b
    var old_r = a

    while r != 0:
        let quotient = old_r / r
        old_r, r = r, old_r - quotient * r
        old_s, s = s, old_s - quotient * s
        old_t, t = t, old_t - quotient * t

    return old_r, old_s, old_t, t, s

def for2(xs, ys, fun): return for xs.length: fun(xs[_], ys[_])

def crt(xs, ys):
    let p = reduce(xs): _a * _b
    var r = 0
    for2(xs,ys) x, y:
        let q = p / x
        let z,s,_t,_qt,_qs = q.extended_gcd(x)
        if z != 1:
            return "ng " + x + " not coprime", 0
        if s < 0: r += y * (s + x) * q
        else:     r += y * s       * q
    return "ok", r % p


def print_crt(xs, ys):
    let msg, res = crt(xs, ys)
    print(msg + " " + res)

print_crt([3,5,7],[2,3,2])
print_crt([11,12,13],[10,4,12])
print_crt([11,22,19],[10,4,9])
print_crt([100,23],[19,0])

Output:

ok 23
ok 1000
ng 11 not coprime 0
ok 1219

Lua

-- Taken from https://www.rosettacode.org/wiki/Sum_and_product_of_an_array#Lua
function prodf(a, ...) return a and a * prodf(...) or 1 end
function prodt(t) return prodf(unpack(t)) end

function mulInv(a, b)
    local b0 = b
    local x0 = 0
    local x1 = 1

    if b == 1 then
        return 1
    end

    while a > 1 do
        local q = math.floor(a / b)
        local amb = math.fmod(a, b)
        a = b
        b = amb
        local xqx = x1 - q * x0
        x1 = x0
        x0 = xqx
    end

    if x1 < 0 then
        x1 = x1 + b0
    end

    return x1
end

function chineseRemainder(n, a)
    local prod = prodt(n)

    local p
    local sm = 0
    for i=1,#n do
        p = prod / n[i]
        sm = sm + a[i] * mulInv(p, n[i]) * p
    end

    return math.fmod(sm, prod)
end

n = {3, 5, 7}
a = {2, 3, 2}
io.write(chineseRemainder(n, a))

Output:

M2000 Interpreter

Translation of: C

Function ChineseRemainder(n(), a()) {
	Function mul_inv(a, b) {
		if b==1 then =1 : exit
		b0=b
		x1=1 : x0=0
		while a>1
			q=a div b
			t=b : b=a mod b: a=t
			t=x0: x0=x1-q*x0:x1=t
		end while
		if x1<0 then x1+=b0
		=x1
	}
	def p, i, prod=1, sum
	for i=0 to  len(n())-1 {prod*=n(i)}
	for i=0 to  len(a())-1
		p=prod div n(i)
		sum+=a(i)*mul_inv(p, n(i))*p
	next
	=sum mod prod
}
Print ChineseRemainder((3,5,7), (2,3,2))

Output:

Maple

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

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

Mathematica / Wolfram Language

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

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

MATLAB / Octave

function f = chineseRemainder(r, m)
  s = prod(m) ./ m;
  [~, t] = gcd(s, m);
  f = s .* t * r';

Output:

>> chineseRemainder([2 3 2], [3 5 7])
 ans = 23

Modula-2

MODULE CRT;
FROM FormatString IMPORT FormatString;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;

PROCEDURE WriteInt(n : INTEGER);
VAR buf : ARRAY[0..15] OF CHAR;
BEGIN
    FormatString("%i", buf, n);
    WriteString(buf)
END WriteInt;

PROCEDURE MulInv(a,b : INTEGER) : INTEGER;
VAR
    b0,x0,x1,q,amb,xqx : INTEGER;
BEGIN
    b0 := b;
    x0 := 0;
    x1 := 1;

    IF b=1 THEN
        RETURN 1
    END;

    WHILE a>1 DO
        q := a DIV b;
        amb := a MOD b;
        a := b;
        b := amb;
        xqx := x1 - q * x0;
        x1 := x0;
        x0 := xqx
    END;

    IF x1<0 THEN
        x1 := x1 + b0
    END;

    RETURN x1
END MulInv;

PROCEDURE ChineseRemainder(n,a : ARRAY OF INTEGER) : INTEGER;
VAR
    i : CARDINAL;
    prod,p,sm : INTEGER;
BEGIN
    prod := n[0];
    FOR i:=1 TO HIGH(n) DO
        prod := prod * n[i]
    END;

    sm := 0;
    FOR i:=0 TO HIGH(n) DO
        p := prod DIV n[i];
        sm := sm + a[i] * MulInv(p, n[i]) * p
    END;

    RETURN sm MOD prod
END ChineseRemainder;

TYPE TA = ARRAY[0..2] OF INTEGER;
VAR n,a : TA;
BEGIN
    n := TA{3, 5, 7};
    a := TA{2, 3, 2};
    WriteInt(ChineseRemainder(n, a));
    WriteLn;

    ReadChar
END CRT.

Output:

Nim

Translation of: C

proc mulInv(a0, b0: int): 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:

OCaml

This is without the Jane Street Ocaml Core Library.

exception Modular_inverse
let inverse_mod a = function
  | 1 -> 1
  | b -> let rec inner a b x0 x1 =
           if a <= 1 then x1
           else if  b = 0 then raise Modular_inverse
           else inner b (a mod b) (x1 - (a / b) * x0) x0 in
         let x = inner a b 0 1 in
         if x < 0 then x + b else x

let chinese_remainder_exn congruences = 
  let mtot = congruences
             |> List.map (fun (_, x) -> x)
             |> List.fold_left ( *) 1 in
   (List.fold_left (fun acc (r, n) ->
                  acc + r * inverse_mod (mtot / n) n * (mtot / n)
                ) 0 congruences)
             mod mtot

let chinese_remainder congruences =
   try Some (chinese_remainder_exn congruences)
   with modular_inverse -> None

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

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:

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]))))
}

Pascal

A console application in Free Pascal, created with the Lazarus IDE.

Follows the Perl examples: (1) expresses each solution as a residue class; (2) if the moduli are not pairwise coprime, still finds a solution if there is one.

// Rosetta Code task "Chinese remainder theorem".
program ChineseRemThm;
uses SysUtils;
type TIntArray = array of integer;

// Defining EXTRA adds optional explanatory code
{$DEFINE EXTRA}

// Return (if possible) a residue res_out that satifies
//   res_out = res1 modulo mod1,  res_out = res2 modulo mod2.
// Return mod_out = LCM( mod1, mod2), or mod_out = 0 if there's no solution.
procedure Solve2( const res1, res2, mod1, mod2 : integer;
                  out res_out, mod_out : integer);
var
  a, c, d, k, m, m1, m2, r, temp : integer;
  p, p_prev : integer;
{$IFDEF EXTRA}
  q, q_prev : integer;
{$ENDIF}
begin
  if (mod1 = 0) or (mod2 = 0) then
    raise SysUtils.Exception.Create( 'Solve2: Modulus cannot be 0');
  m1 := Abs( mod1);
  m2 := Abs( mod2);
  // Extended Euclid's algorithm for HCF( m1, m2), except that only one
  //   of the Bezout coefficients is needed (here p, could have used q)
  c := m1; d := m2;
  p :=0;  p_prev := 1;
{$IFDEF EXTRA}
  q := 1; q_prev := 0;
{$ENDIF}
  a := 0;
  while (d > 0) do begin
    temp := p_prev - a*p;  p_prev := p;  p := temp;
  {$IFDEF EXTRA}
    temp := q_prev - a*q;  q_prev := q;  q := temp;
  {$ENDIF}
    a := c div d;
    temp := c - a*d;  c := d;  d := temp;
  end;
  // Here with c = HCF( m1, m2)
{$IFDEF EXTRA}
  Assert( c = p*m2 + q*m1); // p and q are the Bezout coefficients
{$ENDIF}
  // A soution exists iff c divides (res2 - res1)
  k := (res2 - res1) div c;
  if res2 - res1 <> k*c then begin
    res_out := 0;  mod_out := 0; // indicate that there's no xolution
  end
  else begin
    m := (m1 div c) * m2; // m := LCM( m1, m2)
    r:= res2 - k*p*m2;    // r := a solution modulo m
{$IFDEF EXTRA}
    Assert( r = res1 + k*q*m1); // alternative formula in terms of q
{$ENDIF}
    // Return the solution in the range 0..(m - 1)
    // Don't trust the compiler with a negative argument to mod
    if (r >= 0) then r := r mod m
    else begin
      r := (-r) mod m;
      if (r > 0) then r := m - r;
    end;
    res_out := r;  mod_out := m;
  end;
end;

// Return (if possible) a residue res_out that satifies
//   res_out = res_array[j] modulo mod_array[j], for j = 0..High(res_array).
// Return mod_out = LCM of the moduli, or mod_out = 0 if there's no solution.
procedure SolveMulti( const res_array, mod_array : TIntArray;
                      out res_out, mod_out : integer);
var
  count, k, m, r : integer;
begin
  count := Length( mod_array);
  if count <> Length( res_array) then
    raise SysUtils.Exception.Create( 'Arrays are different sizes')
  else if count = 0 then
    raise SysUtils.Exception.Create( 'Arrays are empty');
  k := 1;
  m := mod_array[0];  r := res_array[0];
  while (k < count) and (m > 0) do begin
    Solve2( r, res_array[k], m, mod_array[k], r, m);
    inc(k);
  end;
  res_out := r;  mod_out := m;
end;

// Cosmetic to turn an integer array into a string for printout.
function ArrayToString( a : TIntArray) : string;
var
  j : integer;
begin
  result := '[';
  for j := 0 to High(a) do begin
    result := result + SysUtils.IntToStr(a[j]);
    if j < High(a) then result := result + ', '
                   else result := result + ']';
  end;
end;

// For the passed-in res_array and mod_array, show the solution
//   found by SolveMulti (above), or state that there's no solution.
procedure ShowSolution( const res_array, mod_array : TIntArray);
var
  mod_out, res_out : integer;
begin
  SolveMulti( res_array, mod_array, res_out, mod_out);
  Write( ArrayToString( res_array) + ' mod '
       + ArrayToString( mod_array) + ' --> ');
  if mod_out = 0 then
    WriteLn( 'No solution')
  else
    WriteLn( SysUtils.Format( '%d mod %d', [res_out, mod_out]));
end;

// Main routine. Examples for Rosetta Code task.
begin
  ShowSolution([2, 3, 2], [3, 5, 7]);
  ShowSolution([3, 5, 7], [2, 3, 2]);
  ShowSolution([10, 4, 12], [11, 12, 13]);
  ShowSolution([1, 2, 3, 4], [5, 7, 9, 11]);
  ShowSolution([11, 22, 19], [10, 4, 9]);
  ShowSolution([2328, 410], [16256, 5418]);
  ShowSolution([19, 0], [100, 23]);
end.

Output:

[2, 3, 2] mod [3, 5, 7] --> 23 mod 105
[3, 5, 7] mod [2, 3, 2] --> 5 mod 6
[10, 4, 12] mod [11, 12, 13] --> 1000 mod 1716
[1, 2, 3, 4] mod [5, 7, 9, 11] --> 1731 mod 3465
[11, 22, 19] mod [10, 4, 9] --> No solution
[2328, 410] mod [16256, 5418] --> 28450328 mod 44037504
[19, 0] mod [100, 23] --> 1219 mod 2300

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

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

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

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

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 lcm/;
say chinese( [2328,16256], [410,5418] ), " mod ", lcm(16256,5418);

Output:

28450328 mod 44037504

Phix

Translation of: C

Uses the function mul_inv() from Modular_inverse#Phix (reproduced below)

function mul_inv(integer a, n)
    if n<0 then n = -n end if
    if a<0 then a = n - mod(-a,n) end if
    integer t = 0,  nt = 1,
            r = n,  nr = a;   
    while nr!=0 do
        integer q = floor(r/nr)
        {t, nt} = {nt, t-q*nt}
        {r, nr} = {nr, r-q*nr}
    end while
    if r>1 then return "a is not invertible" end if
    if t<0 then t += n end if
    return t
end function

function chinese_remainder(sequence n, a)
    integer p, prod = 1, tot = 0;
    for i=1 to length(n) do prod *= n[i] end for
    for i=1 to length(n) do
        p = prod / n[i];
        object m = mul_inv(p, n[i])
        if string(m) then return "fail" end if
        tot += a[i] * m * p;
    end for
    return mod(tot,prod)
end function
 
?chinese_remainder({3,5,7},{2,3,2})
?chinese_remainder({11,12,13},{10,4,12})
?chinese_remainder({11,22,19},{10,4,9})
?chinese_remainder({100,23},{19,0})

Output:

23
1000
"fail"
1219

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)

Prolog

Created with SWI Prolog.

product(A, B, C) :- C is A*B.

pair(X, Y, X-Y).

egcd(_, 0, 1, 0) :- !.
egcd(A, B, X, Y) :-
    divmod(A, B, Q, R),
    egcd(B, R, S, X),
    Y is S - Q*X.

modinv(A, B, X) :-
    egcd(A, B, X, Y),
    A*X + B*Y =:= 1.

crt_fold(A, M, P, R0, R1) :- % system of equations of (x = a) (mod m); p = M/m
    modinv(P, M, Inv),
    R1 is R0 + A*Inv*P.

crt(Pairs, N) :-
    maplist(pair, As, Ms, Pairs),
    foldl(product, Ms, 1, M),
    maplist(divmod(M), Ms, Ps, _), % p(n) <- M/m(n)
    foldl(crt_fold, As, Ms, Ps, 0, N0),
    N is N0 mod M.

Output:

?- crt([2-3, 3-5, 2-7], X).
X = 23.

?- crt([10-11, 4-12, 12-13], X).
X = 1000.

?- crt([2-3, 1-6], X).
false.

gprolog version

Although consult under gprolog didn't complain, I had a problem with "foldl" above. Therefore, I had to write another version of the Chinese Remainder Therorem

/* Chinese remainder Theorem: Input chinrest([2,3,2], [3,5,7], R).  -----> R == 23
                                 or chinrest([2,3], [5,13], R). ---------> R == 42
    Written along the lines of "Introduction to Algorithms" by 
    Thomas Cormen
    Charles Leiserson
    Ronald Rivest
compiled with gprolog 1.4.5 (64 Bits)
*/

chinrest(A, N, X) :- 
    sort(N),                
    prime(N,Nn), !, lenok(A, Nn),                         /* test as to whether the ni are primes */
    product(Nn,P), !,                                     /* P is the product of the ni */
    milist(P, Nn, Mi),                                    /* The Mi List: mi = n/ni */
    cilist(Mi, Nn, Ci),                                   /* The first Ci List: mi-1 mod ni */
    mult_lists(Mi, Ci, Ac),                               /* The ci List :mi*(mi-1 mod ni) */
    mult_lists(Ac, A, Ad),                                /* The ai*ci List */
    sum_list(Ad, S),                                      /* Sum of the ai*cis */
    X is S mod P, ! .                                     /* a is (a1c1 + ... +akck) mod n */

prime([X|Ys], Zs) :- fd_not_prime(X), !, prime(Ys,Zs).    /* sift the primes of [list] */
prime([Y|Ys], [Y|Zs]) :- fd_prime(Y), !, prime(Ys,Zs).
prime([],[]).

product([], 0).                                           /* n1.n2.n3. ... .ni. ... .nk */
product([H|T], P) :- product_1(T, H, P).

product_1([], P, P).
product_1([H|T], H0, P) :- product_1(T, H, P0), P is P0 * H0.
    
lenok(A, N) :- length(A, X), length(N, Y), X=:=Y.
lenok(_, _) :- write('Please enter equal length prime numbers only'), fail.

cilist(Mi, Ni, Ci) :- maplist( (modinv), Mi, Ni, Ci).       /* generate the Cis */

mult_lists(Ai, Ci, Ac) :- maplist( (pro), Ai, Ci, Ac).      /* The mi*ci */
pro(X, Y, Z) :- Z is X * Y.

milist(_, [],[]).
milist(P, [H|T],[X|Y]) :- X is truncate(P/H), milist(P, T, Y).

modinv(A, B, N) :- eeuclid(A, B, P, _, GCD), 
        GCD =:= 1,
        N is P mod B.
        
eeuclid(A,B,P,S,GCD) :-
   A >= B,
   a_b_p_s_(A,B,P,S,1-0,0-1,GCD),
   GCD is A*P + B*S.

eeuclid(A,B,P,S,GCD) :-
   A < B,
   a_b_p_s_(B,A,S,P,1-0,0-1,GCD);
   GCD is A*P + B*S.

a_b_p_s_(A,0,P1,S1,P1-_P2,S1-_S2,A).
a_b_p_s_(A,B,P,S,P1-P2,S1-S2,GCD) :-
   B > 0,
   A > B,
   Q is truncate(A/B),
   B1 is A mod B,
   P3 is P1-(Q*P2),
   S3 is S1-(Q*S2),
   a_b_p_s_(B,B1,P,S,P2-P3,S2-S3,GCD).

Output:

?- chinrest([2,3,2], [3,5,7], X).
X = 23

?- chinrest([2,3], [5,13], X).
X = 42

PureBasic

EnableExplicit
DisableDebugger
DataSection
  LBL_n1:
  Data.i 3,5,7    
  LBL_a1:
  Data.i 2,3,2    
  LBL_n2:
  Data.i 11,12,13
  LBL_a2:
  Data.i 10,4,12
  LBL_n3:
  Data.i 10,4,9
  LBL_a3:
  Data.i 11,22,19
EndDataSection

Procedure ErrorHdl()
  Print(ErrorMessage())
  Input()
EndProcedure

Macro PrintData(n,a)
  Define Idx.i=0
  Print("[")
  While n+SizeOf(Integer)*Idx<a
    Print("( ")
    Print(Str(PeekI(a+SizeOf(Integer)*Idx)))
    Print(" . ")
    Print(Str(PeekI(n+SizeOf(Integer)*Idx)))
    Print(" )")
    Idx+1
  Wend
  Print(~"]\nx = ")
EndMacro

Procedure.i Produkt_n(n_Adr.i,a_Adr.i)
  Define p.i=1
  While n_Adr<a_Adr
    p*PeekI(n_Adr)
    n_Adr+SizeOf(Integer)
  Wend
  ProcedureReturn p
EndProcedure

Procedure.i Eval_x1(a.i,b.i)
  Define b0.i=b, x0.i=0, x1.i=1, q.i, t.i
  If b=1 : ProcedureReturn x1 : EndIf  
  While a>1
    q=Int(a/b)
    t=b : b=a%b : a=t
    t=x0 : x0=x1-q*x0 : x1=t
  Wend
  If x1<0 : ProcedureReturn x1+b0 : EndIf
  ProcedureReturn x1  
EndProcedure

Procedure.i ChineseRem(n_Adr.i,a_Adr.i)
  Define prod.i=Produkt_n(n_Adr,a_Adr), a.i, b.i, p.i, Idx.i=0, sum.i
  While n_Adr+SizeOf(Integer)*Idx<a_Adr  
    b=PeekI(n_Adr+SizeOf(Integer)*Idx)
    p=Int(prod/b) : a=p 
    sum+PeekI(a_Adr+SizeOf(Integer)*Idx)*Eval_x1(a,b)*p
    Idx+1
  Wend
  ProcedureReturn sum%prod
EndProcedure

OnErrorCall(@ErrorHdl())
OpenConsole("Chinese remainder theorem")
PrintData(?LBL_n1,?LBL_a1)
PrintN(Str(ChineseRem(?LBL_n1,?LBL_a1)))
PrintData(?LBL_n2,?LBL_a2)
PrintN(Str(ChineseRem(?LBL_n2,?LBL_a2)))
PrintData(?LBL_n3,?LBL_a3)
PrintN(Str(ChineseRem(?LBL_n3,?LBL_a3)))
Input()

Output:

[( 2 . 3 )( 3 . 5 )( 2 . 7 )]
x = 23
[( 10 . 11 )( 4 . 12 )( 12 . 13 )]
x = 1000
[( 11 . 10 )( 22 . 4 )( 19 . 9 )]
x = Division by zero

Python

Procedural

Python 2.7

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

Python 3.6

# Python 3.6
from functools import reduce
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:

Functional

Using an option type to represent the possibility that there may or may not be a solution for any given pair of input lists.

(Note that the procedural versions above both fail with a ZeroDivisionError on inputs for which no solution is found).

Translation of: Haskell

Works with: Python version 3.7

'''Chinese remainder theorem'''

from operator import (add, mul)
from functools import reduce


# cnRemainder :: [Int] -> [Int] -> Either String Int
def cnRemainder(ms):
    '''Chinese remainder theorem.
       (moduli, residues) -> Either explanation or solution
    '''
    def go(ms, rs):
        mp = numericProduct(ms)
        cms = [(mp // x) for x in ms]

        def possibleSoln(invs):
            return Right(
                sum(map(
                    mul,
                    cms, map(mul, rs, invs)
                )) % mp
            )
        return bindLR(
            zipWithEither(modMultInv)(cms)(ms)
        )(possibleSoln)

    return lambda rs: go(ms, rs)


# modMultInv :: Int -> Int -> Either String Int
def modMultInv(a, b):
    '''Modular multiplicative inverse.'''
    x, y = eGcd(a, b)
    return Right(x) if 1 == (a * x + b * y) else (
        Left('no modular inverse for ' + str(a) + ' and ' + str(b))
    )


# egcd :: Int -> Int -> (Int, Int)
def eGcd(a, b):
    '''Extended greatest common divisor.'''
    def go(a, b):
        if 0 == b:
            return (1, 0)
        else:
            q, r = divmod(a, b)
            (s, t) = go(b, r)
            return (t, s - q * t)
    return go(a, b)


# TEST ----------------------------------------------------
# main :: IO ()
def main():
    '''Tests of soluble and insoluble cases.'''

    print(
        fTable(
            __doc__ + ':\n\n         (moduli, residues) -> ' + (
                'Either solution or explanation\n'
            )
        )(repr)(
            either(compose(quoted("'"))(curry(add)('No solution: ')))(
                compose(quoted(' '))(repr)
            )
        )(uncurry(cnRemainder))([
            ([10, 4, 12], [11, 12, 13]),
            ([11, 12, 13], [10, 4, 12]),
            ([10, 4, 9], [11, 22, 19]),
            ([3, 5, 7], [2, 3, 2]),
            ([2, 3, 2], [3, 5, 7])
        ])
    )


# GENERIC -------------------------------------------------

# Left :: a -> Either a b
def Left(x):
    '''Constructor for an empty Either (option type) value
       with an associated string.'''
    return {'type': 'Either', 'Right': None, 'Left': x}


# Right :: b -> Either a b
def Right(x):
    '''Constructor for a populated Either (option type) value'''
    return {'type': 'Either', 'Left': None, 'Right': x}
# any :: (a -> Bool) -> [a] -> Bool


def any_(p):
    '''True if p(x) holds for at least
       one item in xs.'''
    def go(xs):
        for x in xs:
            if p(x):
                return True
        return False
    return lambda xs: go(xs)


# bindLR (>>=) :: Either a -> (a -> Either b) -> Either b
def bindLR(m):
    '''Either monad injection operator.
       Two computations sequentially composed,
       with any value produced by the first
       passed as an argument to the second.'''
    return lambda mf: (
        mf(m.get('Right')) if None is m.get('Left') else m
    )


# compose (<<<) :: (b -> c) -> (a -> b) -> a -> c
def compose(g):
    '''Right to left function composition.'''
    return lambda f: lambda x: g(f(x))


# curry :: ((a, b) -> c) -> a -> b -> c
def curry(f):
    '''A curried function derived
       from an uncurried function.'''
    return lambda a: lambda b: f(a, b)


# either :: (a -> c) -> (b -> c) -> Either a b -> c
def either(fl):
    '''The application of fl to e if e is a Left value,
       or the application of fr to e if e is a Right value.'''
    return lambda fr: lambda e: fl(e['Left']) if (
        None is e['Right']
    ) else fr(e['Right'])


# fTable :: String -> (a -> String) ->
#                     (b -> String) -> (a -> b) -> [a] -> String
def fTable(s):
    '''Heading -> x display function ->
                 fx display function ->
          f -> value list -> tabular string.'''
    def go(xShow, fxShow, f, xs):
        w = max(map(compose(len)(xShow), xs))
        return s + '\n' + '\n'.join([
            xShow(x).rjust(w, ' ') + (' -> ') + fxShow(f(x))
            for x in xs
        ])
    return lambda xShow: lambda fxShow: lambda f: lambda xs: go(
        xShow, fxShow, f, xs
    )


# numericProduct :: [Num] -> Num
def numericProduct(xs):
    '''The arithmetic product of all numbers in xs.'''
    return reduce(mul, xs, 1)


# partitionEithers :: [Either a b] -> ([a],[b])
def partitionEithers(lrs):
    '''A list of Either values partitioned into a tuple
       of two lists, with all Left elements extracted
       into the first list, and Right elements
       extracted into the second list.
    '''
    def go(a, x):
        ls, rs = a
        r = x.get('Right')
        return (ls + [x.get('Left')], rs) if None is r else (
            ls, rs + [r]
        )
    return reduce(go, lrs, ([], []))


# quoted :: Char -> String -> String
def quoted(c):
    '''A string flanked on both sides
       by a specified quote character.
    '''
    return lambda s: c + s + c


# uncurry :: (a -> b -> c) -> ((a, b) -> c)
def uncurry(f):
    '''A function over a tuple,
       derived from a curried function.'''
    return lambda xy: f(xy[0])(xy[1])


# zipWithEither :: (a -> b -> Either String  c)
#            -> [a] -> [b] -> Either String [c]
def zipWithEither(f):
    '''Either a list of results if f succeeds with every pair
       in the zip of xs and ys, or an explanatory string
       if any application of f returns no result.
    '''
    def go(xs, ys):
        ls, rs = partitionEithers(map(f, xs, ys))
        return Left(ls[0]) if ls else Right(rs)
    return lambda xs: lambda ys: go(xs, ys)


# MAIN ---
if __name__ == '__main__':
    main()

Output:

Chinese remainder theorem:

         (moduli, residues) -> Either solution or explanation

([10, 4, 12], [11, 12, 13]) -> 'No solution: no modular inverse for 48 and 10'
([11, 12, 13], [10, 4, 12]) ->  1000 
 ([10, 4, 9], [11, 22, 19]) -> 'No solution: no modular inverse for 36 and 10'
     ([3, 5, 7], [2, 3, 2]) ->  23 
     ([2, 3, 2], [3, 5, 7]) -> 'No solution: no modular inverse for 6 and 2'

R

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:

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

Output:

Raku

(formerly Perl 6)

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:

REXX

algebraic

/*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 0                                     /*stick a fork in it,  we're all done. */
      end      /*x*/

say 'no solution found.'                         /*oops, announce that solution ¬ found.*/

output when using the default inputs:

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

found a solution with X= 23

congruences sets

/*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             /*build a (array) list of modulo values*/
        end   /*N*/
      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                     /*Is in all sets ? */
    end   /*v*/
  say 'found a solution with X='    x            /*display one possible solution.       */
  exit 0                                         /*stick a fork in it,  we're all done. */
  end   /*x*/

say 'no solution found.'                         /*oops, announce that solution ¬ found.*/

output is identical to the 1^st REXX version.

RPL

Translation of: Python

Works with: Halcyon Calc version 4.2.9

RPL code

Comment

  ≪ 
   DUP ROT 1 R→C ROT 0 R→C
   WHILE DUP RE REPEAT
      OVER RE OVER RE / FLOOR
      OVER * NEG ROT +
   END
   DROP C→R ROT MOD
   SWAP 1 == SWAP 0 IFTE
≫ ‘MODINV’ STO

≪ → n a 
  ≪ 0
     1 1 n SIZE FOR j n j GET * NEXT
     1 n SIZE FOR j
        DUP n j GET / FLOOR
        ROT OVER n j GET MODINV a j GET * ROT * +
        SWAP NEXT
     MOD
≫ ≫ ‘NCHREM’ STO

MODINV ( a b → x )  with ax = 1 mod b
3:b   2:(r,u) ← (a,1)   1:(r',u') ← (b,0)
While r' ≠ 0
    q ← r // r'
   (r - q*r', u - q*u')

Forget (r',u') and calculate u mod b
Test r and return zero if a and b are not co-prime


NCHREM ( n a → remainder )
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
    // reorder stack for next iteration
return sum % prod

RPL 2003 version

Works with: HP version 49

≪ → a n 
  ≪ a 1 GET n 1 GET →V2
     2 a SIZE FOR j
        a j GET n j GET →V2 ICHINREM
        NEXT 
     V→ MOD
≫ ≫ 'NCHREM' STO

{ 2 3 2 } { 3 5 7 } NCHREM

Output:

1: 23.

Ruby

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 = 0, 1
  while remainder != 0
    last_remainder, (quotient, remainder) = remainder, last_remainder.divmod(remainder)
    x, last_x = last_x - quotient*x, x
  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

Rust

fn egcd(a: i64, b: i64) -> (i64, i64, i64) {
    if a == 0 {
        (b, 0, 1)
    } else {
        let (g, x, y) = egcd(b % a, a);
        (g, y - (b / a) * x, x)
    }
}

fn mod_inv(x: i64, n: i64) -> Option<i64> {
    let (g, x, _) = egcd(x, n);
    if g == 1 {
        Some((x % n + n) % n)
    } else {
        None
    }
}

fn chinese_remainder(residues: &[i64], modulii: &[i64]) -> Option<i64> {
    let prod = modulii.iter().product::<i64>();

    let mut sum = 0;

    for (&residue, &modulus) in residues.iter().zip(modulii) {
        let p = prod / modulus;
        sum += residue * mod_inv(p, modulus)? * p
    }

    Some(sum % prod)
}

fn main() {
    let modulii = [3,5,7];
    let residues = [2,3,2];

    match chinese_remainder(&residues, &modulii) {
        Some(sol) => println!("{}", sol),
        None      => println!("modulii not pairwise coprime")
    }

}

Scala

Output:

Best seen running in your browser either by ScalaFiddle (ES aka JavaScript, non JVM) or Scastie (remote JVM).

import scala.util.{Success, Try}

object ChineseRemainderTheorem extends App {

  def chineseRemainder(n: List[Int], a: List[Int]): Option[Int] = {
    require(n.size == a.size)
    val prod = n.product

    def iter(n: List[Int], a: List[Int], sm: Int): Int = {
      def mulInv(a: Int, b: Int): Int = {
        def loop(a: Int, b: Int, x0: Int, x1: Int): Int = {
          if (a > 1) loop(b, a % b, x1 - (a / b) * x0, x0) else x1
        }

        if (b == 1) 1
        else {
          val x1 = loop(a, b, 0, 1)
          if (x1 < 0) x1 + b else x1
        }
      }

      if (n.nonEmpty) {
        val p = prod / n.head

        iter(n.tail, a.tail, sm + a.head * mulInv(p, n.head) * p)
      } else sm
    }

    Try {
      iter(n, a, 0) % prod
    } match {
      case Success(v) => Some(v)
      case _          => None
    }
  }

  println(chineseRemainder(List(3, 5, 7), List(2, 3, 2)))
  println(chineseRemainder(List(11, 12, 13), List(10, 4, 12)))
  println(chineseRemainder(List(11, 22, 19), List(10, 4, 9)))

}

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;

Output:

Sidef

Translation of: Raku

func chinese_remainder(*n) {
    var N = n.prod
    func (*a) {
        n.range.sum { |i|
            var p = (N / n[i])
            a[i] * p.invmod(n[i]) * p
        } % N
    }
}

say chinese_remainder(3, 5, 7)(2, 3, 2)

Output:

SQL

Works with: SQLite version 3.34.0

create temporary table inputs(remainder int, modulus int);

insert into inputs values (2, 3), (3, 5), (2, 7);

with recursive

-- Multiply out the product of moduli
multiplication(idx, product) as (
    select 1, 1

    union all

    select
        multiplication.idx+1,
        multiplication.product * inputs.modulus
    from
        multiplication,
        inputs
    where
        inputs.rowid = multiplication.idx
),

-- Take the final value from the product table
product(final_value) as (
    select max(product) from multiplication
),

-- Calculate the multiplicative inverse from each equation
multiplicative_inverse(id, a, b, x, y) as (
    select
        inputs.modulus,
        product.final_value / inputs.modulus,
        inputs.modulus,
        0,
        1
    from
        inputs,
        product

    union all

    select
        id,
        b, a%b,
        y - (a/b)*x, x
    from
        multiplicative_inverse
    where
        a>0
)
-- Combine residues into final answer
select
    sum(
        (y % inputs.modulus) * inputs.remainder * (product.final_value / inputs.modulus)
    ) % product.final_value
  from
    multiplicative_inverse, product, inputs
  where
    a=1 and multiplicative_inverse.id = inputs.modulus;

Output:

Swift

import Darwin

/*
 * Function: euclid
 * Usage: (r,s) = euclid(m,n)
 * --------------------------
 * The extended Euclidean algorithm subsequently performs
 * Euclidean divisions till the remainder is zero and then
 * returns the Bézout coefficients r and s.
 */

func euclid(_ m:Int, _ n:Int) -> (Int,Int) {
    if m % n == 0 {
        return (0,1)
    } else {
        let rs = euclid(n % m, m)
        let r = rs.1 - rs.0 * (n / m)
        let s = rs.0

        return (r,s)
    }
}

/*
 * Function: gcd
 * Usage: x = gcd(m,n)
 * -------------------
 * The greatest common divisor of two numbers a and b
 * is expressed by ax + by = gcd(a,b) where x and y are
 * the Bézout coefficients as determined by the extended
 * euclidean algorithm.
 */

func gcd(_ m:Int, _ n:Int) -> Int {
    let rs = euclid(m, n)
    return m * rs.0 + n * rs.1
}

/*
 * Function: coprime
 * Usage: truth = coprime(m,n)
 * ---------------------------
 * If two values are coprime, their greatest common
 * divisor is 1.
 */

func coprime(_ m:Int, _ n:Int) -> Bool {
    return gcd(m,n) == 1 ? true : false
}

coprime(14,26)
//coprime(2,4)

/*
 * Function: crt
 * Usage: x = crt(a,n)
 * -------------------
 * The Chinese Remainder Theorem supposes that given the
 * integers n_1...n_k that are pairwise co-prime, then for
 * any sequence of integers a_1...a_k there exists an integer
 * x that solves the system of linear congruences:
 *
 *   x === a_1 (mod n_1)
 *   ...
 *   x === a_k (mod n_k)
 */

func crt(_ a_i:[Int], _ n_i:[Int]) -> Int {
    // There is no identity operator for elements of [Int].
    // The offset of the elements of an enumerated sequence
    // can be used instead, to determine if two elements of the same
    // array are the same.
    let divs = n_i.enumerated()
    
    // Check if elements of n_i are pairwise coprime divs.filter{ $0.0 < n.0 }
    divs.forEach{
        n in divs.filter{ $0.0 < n.0 }.forEach{
            assert(coprime(n.1, $0.1))
        }
    }
    
    // Calculate factor N
    let N = n_i.map{$0}.reduce(1, *)
    
    // Euclidean algorithm determines s_i (and r_i)
    var s:[Int] = []
    
    // Using euclidean algorithm to calculate r_i, s_i
    n_i.forEach{ s += [euclid($0, N / $0).1] }
    
    // Solve for x
    var x = 0
    a_i.enumerated().forEach{
        x += $0.1 * s[$0.0] * N / n_i[$0.0]
    }

    // Return minimal solution
    return x % N
}

let a = [2,3,2]
let n = [3,5,7]

let x = crt(a,n)

print(x)

Output:

Tcl

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:

uBasic/4tH

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

VBA

Uses the function mul_inv() from Modular_inverse#VBA

Translation of: Phix

Private Function chinese_remainder(n As Variant, a As Variant) As Variant
    Dim p As Long, prod As Long, tot As Long
    prod = 1: tot = 0
    For i = 1 To UBound(n)
        prod = prod * n(i)
    Next i
    Dim m As Variant
    For i = 1 To UBound(n)
        p = prod / n(i)
        m = mul_inv(p, n(i))
        If WorksheetFunction.IsText(m) Then
            chinese_remainder = "fail"
            Exit Function
        End If
        tot = tot + a(i) * m * p
    Next i
    chinese_remainder = tot Mod prod
End Function
Public Sub re()
    Debug.Print chinese_remainder([{3,5,7}], [{2,3,2}])
    Debug.Print chinese_remainder([{11,12,13}], [{10,4,12}])
    Debug.Print chinese_remainder([{11,22,19}], [{10,4,9}])
    Debug.Print chinese_remainder([{100,23}], [{19,0}])
End Sub

Output:

Visual Basic .NET

Translation of: C#

Module Module1

    Function ModularMultiplicativeInverse(a As Integer, m As Integer) As Integer
        Dim b = a Mod m
        For x = 1 To m - 1
            If (b * x) Mod m = 1 Then
                Return x
            End If
        Next
        Return 1
    End Function

    Function Solve(n As Integer(), a As Integer()) As Integer
        Dim prod = n.Aggregate(1, Function(i, j) i * j)
        Dim sm = 0
        Dim p As Integer
        For i = 0 To n.Length - 1
            p = prod / n(i)
            sm = sm + a(i) * ModularMultiplicativeInverse(p, n(i)) * p
        Next
        Return sm Mod prod
    End Function

    Sub Main()
        Dim n = {3, 5, 7}
        Dim a = {2, 3, 2}

        Dim result = Solve(n, a)

        Dim counter = 0
        Dim maxCount = n.Length - 1
        While counter <= maxCount
            Console.WriteLine($"{result} = {a(counter)} (mod {n(counter)})")
            counter = counter + 1
        End While
    End Sub

End Module

Output:

23 = 2 (mod 3)
23 = 3 (mod 5)
23 = 2 (mod 7)

Wren

Translation of: C

/* returns x where (a * x) % b == 1 */
var mulInv = Fn.new { |a, b|
    if (b == 1) return 1
    var b0 = b
    var x0 = 0
    var x1 = 1
    while (a > 1) {
        var q = (a/b).floor
        var t = b
        b = a % b
        a = t
        t = x0
        x0 = x1 - q*x0
        x1 = t
    }
    if (x1 < 0) x1 = x1 + b0
    return x1
}

var chineseRemainder = Fn.new { |n, a|
    var prod = n.reduce { |acc, i| acc * i }
    var sum = 0
    for (i in 0...n.count) {
        var p = (prod/n[i]).floor
        sum = sum + a[i]*mulInv.call(p, n[i])*p
    }
    return sum % prod
}

var n = [3, 5, 7]
var a = [2, 3, 2]
System.print(chineseRemainder.call(n, a))

Output:

XPL0

Translation of: C

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

func MulInv(A, B);              \Returns X where rem((A*X) / B) = 1
int  A, B;
int  B0, T, Q;
int  X0, X1;
[B0:= B;  X0:= 0;  X1:= 1;
if B = 1 then return 1;
while A > 1 do
    [Q:= A / B;
    T:= B;  B:= rem(A/B);  A:= T;
    T:= X0;  X0:= X1 - Q*X0;  X1:= T;
    ];
if X1 < 0 then X1:= X1 + B0;
return X1;
];

func ChineseRem(N, A, Len);
int  N, A, Len;
int  P, I, Prod, Sum;
[Prod:= 1;  Sum:= 0;
for I:= 0 to Len-1 do Prod:= Prod*N(I);
for I:= 0 to Len-1 do
   [P:= Prod / N(I);
   Sum:= Sum  +  A(I) * MulInv(P,N(I)) * P;
   ];
return rem(Sum/Prod);
];

int N, A;
[N:= [ 3, 5, 7 ];
 A:= [ 2, 3, 2 ];
 IntOut(0, ChineseRem(N, A, 3));  CrLf(0);
]

Output:

zkl

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

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

@@ Line 3,137: / Line 3,137: @@
 say 'no solution found.'                         /*oops, announce that solution ¬ found.*/</syntaxhighlight>
 {{out|output|text=&nbsp; is identical to the 1<sup>st</sup> REXX version.}} <br><br>
+=={{header|RPL}}==
+{{trans|Python}}
+{{works with|Halcyon Calc|4.2.9}}
+{| class="wikitable" ≪
+! RPL code
+! Comment
+|-
+|
+   ≪
+    DUP ROT 1 R→C ROT 0 R→C
+    '''WHILE''' DUP RE '''REPEAT'''
+       OVER RE OVER RE / FLOOR
+       OVER * NEG ROT +
+    '''END'''
+    DROP C→R ROT MOD
+    SWAP 1 == SWAP 0 IFTE
+ ≫ ‘<span style="color:blue">MODINV</span>’ STO
+ ≪ → n a
+   ≪ 0
+1 n SIZE '''FOR''' j n j GET * '''NEXT'''
+n SIZE '''FOR''' j
+         DUP n j GET / FLOOR
+         ROT OVER n j GET <span style="color:blue">MODINV</span> a j GET * ROT * +
+         SWAP '''NEXT'''
+      MOD
+ ≫ ≫ ‘<span style="color:blue">NCHREM</span>’ STO
+|
+ <span style="color:blue">MODINV</span> ''( a b → x )  with ax = 1 mod b''
+:b   2:(r,u) ← (a,1)   1:(r',u') ← (b,0)
+ While r' ≠ 0
+     q ← r // r'
+    (r - q*r', u - q*u')
+ Forget (r',u') and calculate u mod b
+ Test r and return zero if a and b are not co-prime
+ <span style="color:blue">NCHREM</span> ''( n a → remainder )''
+ 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
+     // reorder stack for next iteration
+ return sum % prod
+|}
+'''RPL 2003 version'''
+{{works with|HP|49}}
+ ≪ → a n
+   ≪ a 1 GET n 1 GET →V2
+a SIZE '''FOR''' j
+         a j GET n j GET →V2 ICHINREM
+         '''NEXT'''
+      V→ MOD
+ ≫ ≫ '<span style="color:blue">NCHREM</span>' STO
+ { 2 3 2 } { 3 5 7 } <span style="color:blue">NCHREM</span>
+{{out}}
+<pre>
+: 23.
+</pre>
 =={{header|Ruby}}==
 Brute-force.
@@ Line 3,180: / Line 3,244: @@
 p chinese_remainder([10,4,9], [11,22,19]) #=> nil
 </syntaxhighlight>
 =={{header|Rust}}==
 <syntaxhighlight lang="rust">fn egcd(a: i64, b: i64) -> (i64, i64, i64) {