不定方程式を解く

例えば,不定方程式の自然数の組み合わせを全てみつけたい.

numpyを使えばシンプルに解けるなとふと思った.

import numpy as np

x = np.arange(100)[:, np.newaxis]
y = np.arange(100)

res = 2*x + 5*y == 40
res.nonzero()

(array([ 0, 5, 10, 15, 20], dtype=int64), array([8, 6, 4, 2, 0], dtype=int64))

np.argwhere(res)[1:-1]

array([[ 5, 6],
[10, 4],
[15, 2]], dtype=int64)

x = np.arange(1, 100)[:, np.newaxis]
y = np.arange(1, 100)

res = 2*x + 5*y == 40
np.argwhere(res) + 1

array([[ 5, 6],
[10, 4],
[15, 2]], dtype=int64)

 
 
 
非常にタイムリーな質問をみつけた.

Sympy: get all integral solutions for linear equation – StackOverflow

Sympyで解く方法は分からないけど.

同様に,numpyのブロードキャストを利用すれば解ける.

求める解が正の整数解であれば,

import numpy as np

x = np.arange(10)[:, np.newaxis, np.newaxis]
y = np.arange(10)[:, np.newaxis]
z = np.arange(10)

res = x*411 + y*295 + z*161 == 3200
np.argwhere(res)

array([[4, 2, 6]], dtype=int64)

[-20, 20)の範囲の整数が条件なら,

x = np.arange(-20, 20)[:, np.newaxis, np.newaxis]
y = np.arange(-20, 20)[:, np.newaxis]
z = np.arange(-20, 20)

res = x*411 + y*295 + z*161 == 3200
np.argwhere(res) - 20

array([[-2, 18, -8],
[ 1, 10, -1],
[ 4, 2, 6],
[ 7, -6, 13]], dtype=int64)
 
 
 
追記:

質問のタイトルが変わっていた.

Get all positive integral solutions for a linear equation – StackOverflow

正の整数解であれば,この場合「(x, y, z) = (4, 2, 6)」しかない.
式をみれば,愚直にループを回した所で,10^3程度なのが分かる.
そして,numpyを使うと結構綺麗に求められる.

そして,面白いコメントが.

This is a linear Diophantine equation, so the module diophantine should be used.

線形ディオファントス方程式を解く為に,DiophantineモジュールというものがSympyにはあるらしい.

Diophantine – SymPy 1.1.1 documentation

試してみる.

from sympy import init_printing
from sympy import symbols
from sympy.solvers.diophantine import diophantine
init_printing()


x, y, z = symbols("x, y, z", integer=True)
diophantine(x*411 + y*295 + z*161 - 3200)

{(t0,2627t0+161t1−19200,−4816t0−295t1+35200)}

線形ディオファントス方程式(ex. ax + by = c)について,
diophantine(ax + by – c)という「eq = 0」の形にする.

from sympy.solvers.diophantine import diop_solve

diop_solve(x*411 + y*295 + z*161 - 3200)

(t0,2627t0+161t1−19200,−4816t0−295t1+35200)

もある.セットとタプルという違いはあるけど,「more level below the highest API」って書いてあるけど,じゃあ何が違うのか,どう使い分けるのかはよく分からない.

diop_solveはdiop_linearを内部的に呼び出して解いているらしいが,直接呼び出す事もできる.

from sympy.solvers.diophantine import diop_linear

diop_linear(x*411 + y*295 + z*161 - 3200)

(t0,2627t0+161t1−19200,−4816t0−295t1+35200)

3元1次方程式の整数解は無数にあるが,正の整数解の場合は「4,2,6」の組み合わせだけ.
記号定義の時に,正の整数と定義すれば,解けるのかな?と思ったけど,

x, y, z = symbols("x, y, z", integer=True, positive=True)
diophantine(x*411 + y*295 + z*161 - 3200)

{(t0,2627t0+161t1−19200,−4816t0−295t1+35200)}

変わらない.

2元n次或いは,n元2次等は簡単に求まりそうだけど,
今回のケースを解けそうなのは無さそうだなあ.
 
 
 

import numpy as np

A = np.array(
    [
        [1, 0, 0],
        [2627, 161, -19200],
        [-4816, -295, 35200]
    ]
)
B = np.array([4, 2, 6])

np.linalg.solve(A, B)

array([ 4., 54., 1.])

%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

t0 = np.linspace(0, 10, 10)
t1 = np.linspace(0, 60, 10)
x = t0
y = 2627*t0 + 161*t1 - 19200
z = -4816*t0 - 295*t1 + 35200

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)

FireShot Capture 259 - JupyterLab Alpha Preview - http___localhost_8888_lab
 
 
 

%%timeit
x = np.arange(10)[:, np.newaxis, np.newaxis]
y = np.arange(10)[:, np.newaxis]
z = np.arange(10)

res = x*411 + y*295 + z*161 == 3200
np.argwhere(res)

56.2 µs ± 2.51 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

def func_sol():
    for x in range(1, 10):
        for y in range(1, 10):
            for z in range(1, 10):
                if 411*x+295*y+161*z == 3200:
                    return x, y, z

%timeit func_sol()

110 µs ± 1.27 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

広告
カテゴリー: 未分類 パーマリンク

不定方程式を解く への3件のフィードバック

  1. ピンバック: 整数配列からa + b + c = 0を満たす様な組み合わせを求める | 粉末@それは風のように (日記)

  2. ピンバック: Sympyでディオファントス方程式を解く | 粉末@それは風のように (日記)

  3. ピンバック: 制約最適化問題の実現可能な領域を視覚化したい | 粉末@それは風のように (日記)

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中