中国剰余定理の素朴な実装

南北朝時代に書かれた算術書<b>『孫子算経』</b>には次のような問題があったそうです：
> 3で割ると2余り、5で割ると3余り、7で割ると2余る数は何か

これをみたす最小の自然数は$23$です。

他にも$23,128,233,\cdots$と無限に存在するのですが、実はどれも$105$で割った余りが$23$の整数になっています。さらに$105$という数は、割る数の積になっています：$105=3\times5\times7$。

この問題はどんな数で設定しても解が存在するとは限りません。例えば、「3で割ると1余り、3で割ると2余る数」は明らかに存在しませんし、「3で割ると1余り、6で割ると2余る数」も存在しません。割る数には「どの2数も<b>互いに素</b>（最大公約数が$1$）」という条件が必要なのです。


このような問題を一般の整数に拡張した<font color="red"><b>中国剰余定理</b></font>と呼ばれる定理があります。この記事では中国剰余定理の解説と素朴な実装（高速化等はしていない）を紹介します。


<h2>合同式について</h2>
$n$を$m$で割った余りが$k$であることを

```math
n\equiv k\pmod m
```
と表すことにします（合同式といいます）。例えば、$15$を$7$で割った余りは$1$なので$15\equiv1\pmod7$と表せます。$n\equiv0\pmod 6$は$n$を$6$で割った余りが$0$である（$n$は$6$の倍数である）ことを表します。合同式はイコールと同じように、両辺に同じ数を足したり引いたり掛けたりしても成り立つという性質を持っています（割り算はちょっとした条件が必要）。

また、この記事では以下の事実を既知とします：
$m$と$n$が互いに素$\iff xn\equiv 1\pmod m$ を満たす$x$が存在する。

例えば$m=4$，$n=7$のとき、$x=3$とすれば$3\times 7=21\equiv1\pmod4$となります。
逆に互いに素でない場合、例えば$m=4$，$n=6$とすると、$x$が何であろうと$6x$は偶数であるから、$4$で割った余りの候補は$0$か$2$です。このため$6x\not\equiv1\pmod4$となってしまいます。

<h1>中国剰余定理と証明</h1>

それでは中国剰余定理の紹介と証明です。

<b>定理</b>　$n$個の自然数$m_1, \ldots,m_n$はどの2数も互いに素であるとし、$m=m_1\times\cdots\times m_n$とする。このとき、任意の整数$r_1, \ldots,r_n$に対して

```math
\begin{align*}
a&\equiv r_1\pmod{m_1},\\
a&\equiv r_2\pmod{m_2},\\
 &\vdots\\
a&\equiv r_n\pmod{m_n}\\
\end{align*}
```
をみたす数$a$が$\mod m$で一意に存在する。


<b>証明</b>

各$i=1,\cdots,n$に対して、$m_i^* :=\dfrac{m}{m_i}$とおく。このとき、$m_i^* $は整数であり、$m_i$と$m_i^* $は互いに素である。よって$x_i m_i^* \equiv 1 \pmod{m_i}$を満たす整数$x_i$が存在するので、$r_i x_i m_i^* \equiv r_i \pmod{m_i}$となる。

ここで

```math
a := r_1 x_1 m_1^* + r_2 x_2 m_2^* + \cdots + r_n x_n m_n^*
```

とおく。$j\neq i$であれば$m_j^* \equiv0\pmod{m_i}$なので、任意の$i=1,\cdots,n$に対して

```math
a \equiv 0+\cdots + r_i x_i m_i^* + \cdots + 0 \equiv r_i \pmod{m_i}
```

が成り立つので、$a$が求める解となる。

また、$a$、$a'$が整数解であるとすると

```math
\begin{align*}
a&\equiv a'\pmod{m_1},\\
a&\equiv a'\pmod{m_2},\\
 &\vdots\\
a&\equiv a'\pmod{m_n}\\
\end{align*}
```

であるから

```math
\begin{align*}
a-a'&\equiv 0\pmod{m_1},\\
a-a'&\equiv 0\pmod{m_2},\\
 &\vdots\\
a-a'&\equiv 0\pmod{m_n}\\
\end{align*}
```
となる（つまり$a-a'$は$m_i$の倍数）。$m_i$たちはどの2数も互いに素であったから、$a-a'$は$m$の倍数であることが分かる。これは$a\equiv a'\pmod m$であることを示している。■

$a$の構成がうまいですね。$m_i$で割ったときに1つだけ$r_i$余り、それ以外は余りが$0$となるような数をうまく作っています。


<h1>中国剰余定理の実装</h1>
それでは今の証明の流れを意識した中国剰余定理のコードの紹介です。

```python
import numpy as np
import sys

#長さnのリストからr個取り出したときの組み合わせを取得する関数
def combinations(iterable, r):
    pool = tuple(iterable)
    n = len(pool)
    if r > n:
        return
    indices = list(range(r))
    yield tuple(pool[i] for i in indices)
    while True:
        for i in reversed(range(r)):
            if indices[i] != i + n - r:
                break
        else:
            return
        indices[i] += 1
        for j in range(i+1, r):
            indices[j] = indices[j-1] + 1
        yield tuple(pool[i] for i in indices)


#中国剰余定理（Chinese Remainder Theorem）の実装
def CRT(mod_list, remainder_list):

    m = np.prod(mod_list)
    n = len(mod_list)

    #mod_listの中から互いに素でない自然数の組を探し、あったらプログラム終了
    for pair in combinations(mod_list, 2):
        if np.gcd(pair[0],pair[1])!=1:
            print(f"Error:{pair}の最大公約数が{np.gcd(pair[0],pair[1])}です！")
            sys.exit()

    #証明中に作ったm_i^*のリスト
    mod_list2 = [m//k for k in mod_list]

    ri_xi_mi2_list = []
    for i in range(n):
        mi = mod_list[i]
        #ri（余り）をmi（割る数）よりも大きく設定してしまった場合のために、riをmiで割った余りにしておく
        ri = remainder_list[i]%mi
        mi2 = mod_list2[i]
        #xi*mi2≡1 (mod mi)となるxiを0から探す
        xi = 0
        while (xi*mi2)%mi != 1:
            xi += 1
        ri_xi_mi2_list.append(ri*xi*mi2)

    #足し合わせてmで割った余りを計算
    a = np.sum(ri_xi_mi2_list)%m
    return print(f"{a} (mod {m})")


#割る数リスト（どの2数も互いに素であるように自然数を入力してください）
mod_list = [3, 5, 7]
#余りリスト（リストの長さはmod_listと同じにしてください）
remainder_list = [2, 3, 2]

CRT(mod_list, remainder_list)
```

結果

```python
23 (mod 105)
```

<h1>参考文献など</h1>
- <a  href="https://ja.wikipedia.org/wiki/%E4%B8%AD%E5%9B%BD%E3%81%AE%E5%89%B0%E4%BD%99%E5%AE%9A%E7%90%86" target="_blank" rel="noopener noreferrer">中国の剰余定理 - Wikipedia</a>
- <a  href="https://www.amazon.co.jp/%E6%95%B0%E8%AB%96%E5%BA%8F%E8%AA%AC-%E5%B0%8F%E9%87%8E-%E5%AD%9D/dp/4785310502" target="_blank" rel="noopener noreferrer">数論序説</a>
- <a  href="https://docs.python.org/ja/3/library/itertools.html" target="_blank" rel="noopener noreferrer">itertools --- 効率的なループ実行のためのイテレータ生成関数</a>