Pythonで2元1次不定方程式（ax+by=c）を解く

必要知識
* 高校数学（数A）
* $2 \times 2$行列の積、逆行列、単位行列
* Python 3、NumPy少し
* $a,b$ の最大公約数（greatest common divisor）を $\mathrm{gcd}(a,b)$ と表す。

### 二元一次不定方程式とは
$x,y$ の2つを求めなければならないのに1次式が1つしかないので解は定まらないよ〜ってやつ。

___例題___
二元一次不定方程式 $17 x + 12 y = 1$ の整数解を求めよ。

___解法___
$17 = 12 \cdot 1 + 5，12 = 5 \cdot 2 + 2，5 = 2 \cdot 2 + 1$ だから，

```math:
\begin{align}
  1 &= 5 - 2 \cdot 2 \\
  &= 5 - (12 - 5 \cdot 2) \cdot 2 \\
  &= 12 \cdot (-2) + 5 \cdot 5 \\
  &= 12 \cdot (-2) + (17 - 12 \cdot 1) \cdot 5 \\
  &= 17 \cdot 5 + 12 \cdot (-7)
\end{align}
```

$17x+12y=1$ から $17 \cdot 5 + 12 \cdot (-7)=1$ の辺々を引くと，

```math:
\begin{align}
  17(x-5)+12(y+7)=0
\end{align}
```

$17$ と $12$ は互いに素より，$x-5=12k\ (k\in \mathbb{Z})$ と表せて，

```math:
\begin{align}
  x=5+12k，y=-7-17k \quad (k\in \mathbb{Z})
\end{align}
```

### Pythonで ax+by=1 を解く

#### 方針
1. $a,b$ が互いに素でないなら解なし。互いに素なら解が存在する。[^1]
1. ユークリッドの互除法を遡り、特殊解 $(x_0,y_0)$ を見つける。
1. $a(x-x_0)+b(y-y_0)=0$ より、$x-x_0=bk\ (k\in \mathbb{Z})$ と表すことができる。
1. 解は $(x,y)=(x_0+bk, y_0-ak)$


#### Step1. 最大公約数を求める
$\mathrm{gcd}(a,b)=1$ は互いに素と同値。ユークリッドの互除法を用いる。
$(a,b) \rightarrow (b,r) \rightarrow \cdots  \rightarrow (g,0)$ になるまで繰り返すと、$\mathrm{gcd}(a,b)=g$

```python:
def gcd(a,b):
    while b!=0:
        a,b = b,a%b
    return a
```

a < b であっても1手目が$(a,b) \rightarrow (b,a)$となるだけなので問題ない。
こんな計算しなくても、mathモジュールの`math.gcd(a,b)`やNumPyモジュールの`numpy.gcd(a,b)`を使えばよい。


#### Step2. 特殊解を求める（今回のメイン）

##### 互除法の逆
$ax_0+by_0=1$ について，$a=bq_0+r_0$ を代入すると，

```math:
\begin{align}
  (bq_0+r_0)x_0 + by_0 &= 1 \\
  b(q_0x_0+y_0) + rx_0 &= 1
\end{align}
```

$x_1=q_0x_0+y_0$，$y_1=x_0$ と考えれば，$bx_1+r_0y_1=1$
これを繰り返せば $1x_k+0y_k=1$ に到達し，$(x_k,y_k)=(1,0)$

```math:
\begin{align}
  \begin{pmatrix}
    q_0 & 1 \\
    1 & 0
  \end{pmatrix}
  \begin{pmatrix}
    x_0 \\ y_0
  \end{pmatrix}
  &=
  \begin{pmatrix}
    x_1 \\ y_1
  \end{pmatrix}\\
  \begin{pmatrix}
    x_0 \\ y_0
  \end{pmatrix}
  &=
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_0
  \end{pmatrix}
  \begin{pmatrix}
    x_1 \\ y_1
  \end{pmatrix}\\
  \begin{pmatrix}
    x_0 \\ y_0
  \end{pmatrix}
  &=
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_0
  \end{pmatrix}
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_1
  \end{pmatrix}
  \cdots 
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_{k-1}
  \end{pmatrix}
  \begin{pmatrix}
    x_k \\ y_k
  \end{pmatrix}
\end{align}
```

つまり，

```math:
\begin{align}
A_k=
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_0
  \end{pmatrix}
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_1
  \end{pmatrix}
  \cdots 
  \begin{pmatrix}
    0 & 1 \\ 1 & -q_{k-1}
  \end{pmatrix}
\end{align}
```

を求めれば特殊解が求まる。

##### 特殊解を求めるアルゴリズム方針
入力：$(a,b)$
$A_0\ を\ 2\times 2$ の単位行列とする。

$ax_0+by_0=1$ について，割り算 $a=bq+r$ を用いて $bx_1+ry_1=1$ に変形．

```math:
\begin{align}
A_1=A_0
  \begin{pmatrix}
    0 & 1 \\ 1 & -q
  \end{pmatrix}
  \quad 
\end{align}
```

$(a,b,A_0)$ を $(b,r,A_1)$ に更新する．以下繰り返し．

```math:
\vdots
```

$ax_{k-1}+by_{k-1}=1$ について，割り算 $a=1q+0$ を用いて $1x_k+0y_k=1$ に変形．

```math:
\begin{align}
A_k=A_{k-1}
\begin{pmatrix}
0 & 1 \\ 1 & -q
\end{pmatrix}
\end{align}
```

$(a,b,A_0)$ を $(1,0,A_k)$ に更新したところで終了．

出力は，
```math:
\begin{align}
  \begin{pmatrix}
    x_0 \\ y_0
  \end{pmatrix}=A_k
  \begin{pmatrix}
    1 \\ 0
  \end{pmatrix}=
  \begin{pmatrix}
    Aの11成分 \\ Aの21成分
  \end{pmatrix}
\end{align}
```

※ $a<b$の場合でも，同じ操作で$q=0$となり，$(a,b,A_0) \leftarrow (b,a,A_1)$と入れ替わるだけ。

##### 特殊解をPythonで求める

```python
import numpy as np
def ex_euc(a,b):
    A = np.identity(2, dtype=int) # 2×2の単位行列
    while b!=0: # b==0で終了
        q = a // b
        r = a % b
        A = A @ np.array([[0, 1],
                          [1, -q]])
        a,b = b,r
    return A[0][0], A[1][0]
```
NumPyを使いたくない場合、2×2行列の積だけ追加で定義すればできる。

#### Step3. 解を表示する

```python
def solve_1(a,b):
    if np.gcd(a,b)!=1:
        return "解なし"
    else:
        x_0,y_0 = ex_euc(a,b)
        return f"(x,y)=({x_0}+{b}k, {y_0}-{a}k)"

print(solve_1(a,b))
```


### Pythonで ax+by=c を解く（完成）
__方針__
1. $c$ は $\mathrm{gcd}(a,b)=g$ の倍数でないなら解なし。倍数なら解が存在する。[^2]
1. 両辺を $g$ で割って、$a ^ \prime x+b ^ \prime y=c ^ \prime $ にすると $
a ^ \prime , b ^ \prime$ は互いに素。
1. $a ^ \prime x+b ^ \prime y=1 $ は解をもち、$(x,y)=(x_0+b ^ \prime k, y_0-a ^ \prime k)$
1. $a(c ^ \prime x)+b(c ^ \prime y)=c ^ \prime$ は解をもち、$(c ^ \prime x,c ^ \prime y)=(c ^ \prime x_0+c ^ \prime b ^ \prime k, c ^ \prime y_0-c ^ \prime a ^ \prime k)$

```python
import numpy as np

def ex_euc(a,b):
    A = np.identity(2, dtype=int) # 2×2の単位行列
    while b!=0: # b==0で終了
        q = a // b
        r = a % b
        A = A @ np.array([[0, 1],[1, -q]])
        a,b = b,r
    return A[0][0], A[1][0]

def solve_2(a,b,c):
    g = np.gcd(a,b)
    if c%g!=0:
        return "解なし"
    else:
        a = a // g #a'のこと
        b = b // g #b'のこと
        c = c // g #c'のこと
        x_0,y_0 = ex_euc(a,b)
        return f"(x,y)=({c*x_0}+{c*b}k, {c*y_0}-{c*a}k)"
```
実行してみた。
```python
(a,b,c) = (2,4,2)
print(f"{a}x+{b}y={c}")
print(solve_2(a,b,c))
```
>2x+4y=2
>(x,y)=(1+2k, 0-1k)

0を非表示にした方が見た目いいかも。

```python
(a,b,c) = (3,6,7)
print(f"{a}x+{b}y={c}")
print(solve_2(a,b,c))
```
>3x+6y=7
>解なし

左辺だけ3の倍数なので自明。

```python
(a,b,c) = (981254,652312,4568)
print(f"{a}x+{b}y={c}")
print(solve_2(a,b,c))
```
>981254x+652312y=4568
>(x,y)=(-138239100+744940304k, 207949064-1120592068k)

……本当？


[^1]:「$ax+by=1$ が整数解を持つ $\iff a,b$ は互いに素」
$(\implies)$ $a,b$ が互いに素でないなら $ax+by=(\mathrm{gcd}の倍数)\neq 1$ で整数解なし．
$(\impliedby)$ 互いに素なら，$a,2a,⋯,(b−1)a$ を $b$ で割った余りは全て異なり $1〜(b-1)$ なので，余りが $1$ となるようなものが存在し，$na=bq+1$ から整数解を作れる．
（参考）[二元一次不定方程式と鳩ノ巣原理 - Qiita](https://qiita.com/rime_math/items/954a8771fc6b37021daa)

[^2]:「$ax+by=c$ が整数解を持つ $\iff c$ は $\mathrm{gcd}(a,b)$ の倍数」
$(\implies)$ $c$ が $\mathrm{gcd}(a,b)$ の倍数でないとき，整数解なし．
$(\impliedby)$ $c$ が $\mathrm{gcd}(a,b)=g$ の倍数のとき，$a=a ^ \prime g，b=b ^ \prime g\ (a ^ \prime,b ^ \prime は互いに素)$ と表せて，$a ^ \prime x+b ^ \prime y=1$ が整数解を持つ．これの両辺に $c=kg$ をかけると $a(kx)+b(ky)=c$ ．