Comparison Zig, Julia, Elixir, Javascript on Chinese Remainder Theorem

We have all learnt this theorem one day. A refresher maybe? [wiki](https://en.wikipedia.org/wiki/Chinese_remainder_theorem) or maybe you prefer [this short Standford one](https://crypto.stanford.edu/pbc/notes/numbertheory/crt.html).

You have the [Rosetta code](https://rosettacode.org/wiki/Modular_inverse#Elixir) that explains algorithms.

It is about finding a solution(s) to multiple congruences.
It is a famous result extensively used in Algebra, cryptography and computer science

Why is this interesting? To compute a solution, we must find a modular inverse, so use some new primitives.

We can do:
- brut force. We compute `b x i` for each integer `i` in the range `(0..m-1)` and find an inverse when `rem(b * i,m) == 1`.
> It is one of the implementations in Elixir.
- the famous extended euclidean algorithm. We implement this in the "invMod" function.
> To find the inverse of `a mod m`, you search for integral Bezout coefficients `u,v` such that
`gcd(a,m)=ua +vm`
so that when `a` and `m` are coprime, the gcd is 1, so you have an inverse of `a` modulo `m`.
It is implemented with recursion with Elixir and a "while" loop" in Javascript.
- [Julia](https://docs.julialang.org/en/v1/#man-introduction) implements another algorithm.

For example, the inverse of 3 modulo 7 is 5 because 3x5=15 = 2x7 + 1 = 1 mod 7, so we have 3 x 5 = 1 mod 7 (and the integers modulo 7 form the field Z/7Z, so success was guaranteed).

We will check that the solution of following congruences:
x = 2 mod 3
x = 3 mod 5
x = 2 mod 7

is 23 mod 105 (=3x5x7)

Indeed, 23 = 7x3+2 = 4x5+3 = 3x7+2, but also 105+23=128 is also 2[3] and 3[5] and 2[7] etc

And the theorem states that is it unique modulo 3x5x7 (because the moduli are pairwise coprime...).

```
crt([2, 3, 2], [3, 5, 7])

23
```

## Zig implementation

Learning points:
- If you pass arrays in Zig, you have to pass a reference. Since I use a slice, it is already a pointer with a length (slices have defined length).
- Then variables are immutable by default. If you mutate them, create a new variable.
- internal variables in loops are instantiated on every step (thus you can use `const`).

Run this with `zig build crt.zig`, or build it `zig build-exe ctr.zig` and `./crt`.

```zig
const std = @import("std");
const print = std.debug.print;

fn crt(a: []const i64, m: []const i64) !i64 {
// Compute M = product of all moduli
var M: i64 = 1;
for (m) |mi| {
M *= mi;
}

var result: i64 = 0;
for (a, 0..) |ai, i| {
const Mi = @divFloor(M, m[i]);
const inv = try invMod(Mi, m[i]);
result += ai * Mi * inv;
}

return @rem(result, M);
}

fn invMod(a: i64, m: i64) !i64 {
var av = a;
var mv: i64 = m;
var x0: i64 = 0;
var x1: i64 = 1;

if (m == 1) return 1;

while (av > 1) {
const t: i64 = mv;
const q: i64 = @divFloor(av, mv);
mv = @rem(av, mv);
av = t;

const temp: i64 = x1;
x1 = x0;
x0 = temp - q * x0;
}

x1 = @rem(x1 + m, m);
if (av != 1) return error.NoModInverse; // No modular inverse if gcd != 1
return x1;
}

pub fn main() void {
const a = [_]i64{ 2, 3, 2 };
const m = [_]i64{ 3, 5, 7 };
const result = crt(&a, &m);
if (result) |r| {
print("{}\n", .{r});
} else |err| {
print("{}, No solution\n", .{err});
}
}

```

## Julia implementation

The [Julia](https://juliaacademy.com/courses/?preview=logged_out) code is of course very concise, as it is made for vector computations and parallelism.

In [Julia](https://docs.julialang.org/en/v1/manual/mathematical-operations/#Arithmetic-Operators), the integral division uses the operator `\div tab` : ÷.
The "inverse modulo" is implemented with `invmod`. The algorithm is described [here](https://arxiv.org/pdf/2204.04342).

This can be run with `julia crt.jl` (if you have Julia installed, and saved this file as "crt.jl"!).

Julia

```julia
function crt(a::Array, n::Array)
p = prod(n)
inv = sum(ai * invmod(p ÷ ni, ni) * (p ÷ ni) for (ni, ai) in zip(n, a))
mod(inv, p)
end

@show crt([2, 3, 2], [3, 5, 7])
```

## Elixir implementation

With `Elixir`, we compute the inverse modulo by "brut force" and a recursive version of the "extended euclidean".
A fancy usage of `zip_reduce`, and using `Kernel.div` and `Kernel.rem` (not easy to find in the docs... as they are in the ...Kernel guards!

Run it with `elixir crt.exs`.

Elixir

```elixir
defmodule Ex do
def crt(a, m) do
prod = Enum.product(m)

Enum.zip_reduce([a, m], 0, fn [a_i, m_i], acc ->
b_i = div(prod, m_i)
inv_i = inv_mod_bz(b_i, m_i)
acc + a_i * b_i * inv_i
end)
|> rem(prod)

# extended euclidean (Bezout coefficients)
def inv_mod_bz(p_i, 0), do: {p_i, 1, 0}

def inv_mod_bz(p_i, m_i) do
{d, u, v} = inv_mod_bz(m_i, rem(p_i, m_i))
{d, v, u - div(p_i, m_i) * v}
end

# brut force method
def inv_mod_bf(_, 1), do: 1

def inv_mod_bf(b, m_i) do
# brutforce to compute the modular inverse
inverse =
Enum.reduce(1..m_i, nil, fn j, acc ->
if rem(b * j, m_i) == 1 do
j
else
acc
end
end)

if inverse == nil do
nil
else
rem(inverse, m_i)
end
end
end

Ex.crt([2, 3, 2], [3, 5, 7]) |> dbg()
```

## Javascript implementation

With `Javascript`, we also use the extended euclidean to compute the modular inverse.

Run it with `bun crt.js`

Javascript

```js
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 pi = Math.floor(prod / ni);
sum += ri * pi * mulInv(pi, ni);
}
return sum % prod;
}

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

if (ni === 1) {
return 1;
}

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

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

```

## Python implementation

Recall that integral division uses `//`.. and modulo is simply `%`.

Run it with `python3 crt.py`.

Python

```python
def crt(a,n):
prod = 1
for ni in n:
prod *= ni

sum = 0
for n_i, a_i in zip(n, a):
p = prod // n_i
sum += a_i * inv_mod(p, n_i) * p
return sum % prod

def inv_mod(p, n_i):
b0 = n_i
x0, x1 = 0, 1
if n_i == 1: return 1
while p > 1:
q = p // n_i
p, n_i = n_i, p % n_i
x0, x1 = x1 - q * x0, x0
if x1 < 0: x1 += b0
return x1

print(crt([2, 3, 2], [3, 5, 7]))
```