LoginSignup
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@dolce_itf

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

Posted at

必要知識

  • 高校数学（数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$ だから，

\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$ の辺々を引くと，

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

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

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

Pythonで ax+by=1 を解く

方針

  1. $a,b$ が互いに素でないなら解なし。互いに素なら解が存在する。1
  2. ユークリッドの互除法を遡り、特殊解 $(x_0,y_0)$ を見つける。
  3. $a(x-x_0)+b(y-y_0)=0$ より、$x-x_0=bk\ (k\in \mathbb{Z})$ と表すことができる。
  4. 解は $(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$

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$ を代入すると，

\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)$

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

つまり，

\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$ に変形．

\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)$ に更新する．以下繰り返し．

\vdots

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

\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)$ に更新したところで終了．

出力は，

\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で求める
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. 解を表示する

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
  2. 両辺を $g$ で割って、$a ^ \prime x+b ^ \prime y=c ^ \prime $ にすると $
    a ^ \prime , b ^ \prime$ は互いに素。
  3. $a ^ \prime x+b ^ \prime y=1 $ は解をもち、$(x,y)=(x_0+b ^ \prime k, y_0-a ^ \prime k)$
  4. $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)$
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)"

実行してみた。

(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を非表示にした方が見た目いいかも。

(a,b,c) = (3,6,7)
print(f"{a}x+{b}y={c}")
print(solve_2(a,b,c))

3x+6y=7
解なし

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

(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 ↩

  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$ ． ↩

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?