ベズー方程式の解の公式 その1（ax+by=c）
過去ブログの転載です。
これベズーの等式という名前がついているそうですね。
この形の不定方程式の自然数解(x, y)を求める問題はよくありますが、これを解くのって意外とめんどくさいことが多いのですよね。そこで、2次方程式みたいに解の公式を作ればラクになるんじゃないかと思って作ってみました。結局、とても複雑な式になってしまったのですが……少なくとも覚える気はしませんね。
以後、a, bは互いに素で正の数であるとします。ax－by＝cの形もベズーの等式ですけど、次回にします。
ax＋by＝1の特殊解(X, Y)
ただし、f(a, b)は次式で定義する再帰関数とする。
f(a, b)はax＋by＝1の特殊解Xを与える関数として定義されたものです。ちょっとmodの前にスペースができちゃって紛らわしいのですが、分子のはf(b, a mod b)です。一つ目の変数はb mod bじゃなくてただのbです。
ax＋by＝cの自然数解
nは整数で、a, bが正の数のとき、自然数解は有限個になります。自然数解ではなく単なる整数解の場合は、nの条件を取っ払って整数を入れれば、そのままそれが整数解になります。整数解は無数にあります。
上のnの不等式の上限が－1になる場合、解が存在しないことを意味します。解の個数はこの上限値に1を足したものであると言えます。
結局何がしたかったのかというと、単に解の公式を作りたかっただけです。この式によりすべての自然数解を機械的に導くことが可能ですが、簡単な問題であってもだるい手続きが必要なので、公式の適用はただただだるいです。ただし、機械的に導けるということは、プログラミングするときは有用なんでしょうね。上式を用いれば、簡単なプログラムですべての自然数解を列挙できるようになります。
この解の公式は、a, bが互いに素でなければ機能しません。互いに素になってない場合は両辺を何かで割ってください。割れない場合はどうがんばっても右辺の形にならないので、解なしということになります。
例: 152x＋41y＝30000
152x＋41y＝30000の自然数解を求めましょう。ふつうの手計算で求められるかな？
特殊解を見つけて、整数nを使って一般解を作り、x, yが自然数だからx＞0, y＞0としてnの範囲を絞る流れが一般的です（よね……？）。
それでは、解の公式に代入して解を得てみましょう。この場合は(a, b, c)＝(152, 41, 30000)ですね。再帰関数f(a, b)を求めて、特殊解(X, Y)を求めて、nの範囲を求めて、一般解(x, y)を求めるという手順です。
こんなんやるもんじゃないですね。
再帰関数はf(a, b)のbが1になるとf(a, b)＝0になりますから、そうなるまで再帰し続けるため、このようにすごい入れ子構造になってしまいます。やっていると何となく分かると思いますが、ユークリッドの互除法をただ数式にしただけだったりします。
f(152, 41)＝17でした。コンピュータはこういうの強そうで何よりです。
f(152, 41)を用いて、特殊解(X, Y)がただちに求まります。
あとはこれを用いて、一般解(x, y)が解の公式より得られます。
まずはnの上限値を求めておきます。
0≦n≦4ということなので、解は5個あることになりますね。
これが解ということになります。
実際にn＝0, 1, 2, 3, 4を順番に代入することによって、具体的な解をすべて得られます。
ここまで求めておいてなんですが、こりゃ使い物になりませんね。
やはり手計算に用いるのでなく、定式化用だとか、コンピュータ用にするのが順当でしょう。
導出
まず、ax＋by＝1の特殊解を求める手続きを考えます。以後、a, bは正の数で、gcd(a, b)＝1とします。ユークリッドの互除法をすることにより、a, bの係数を小さくしていくことができて、最終的には1と0になるという性質（ラメの定理というそうです）を用います。不思議なことに、a, bのどっちが大きくても適用できるのですが、イメージしやすいように、ここではa＞bであると考えてください。
a÷b＝qあまりrであったとすると、a＝bq＋rが成り立ちます（これは除法の原理）。すると、a－bq＝rが成り立つことになりますね。ユークリッドの互除法では次に何をするかというと、(a, b)のペアから(b, r)のペアに浮気するわけですね。これを繰り返すとどんどんペアが小さくなります。ここで、小さい方のペア(b, r)からなる不定方程式bx＋ry＝1の解が分かっていたとします。この解をそれぞれx＝X(b, r), y＝Y(b, r)という二つの関数でそれぞれ表わしておきます。すると、当然ですけれどbX(b, r)＋rY(b, r)＝1が成り立つわけですね。
a－bq＝r … ①
bX(b, r)＋rY(b, r)＝1 … ②
この二つの式を見比べて、①の右辺がrであることから、②のrに代入できます。
bX(b, r)＋(a－bq)Y(b, r)＝1
bX(b, r)＋aY(b, r)－bqY(b, r)＝1
aY(b, r)＋b{X(b, r)－qY(b, r)}＝1
これをよく見てみると、ax＋by＝1としたとき
x＝Y(b, r), y＝X(b, r)－qY(b, r)になってるわけですよね。
つまり、小さいペア(b, r)からなる方程式bx＋ry＝1の特殊解が分かっていれば、それを用いて大きいペア(a, b)からなる方程式ax＋by＝1の特殊解も分かるということになります。
なんかかっこつけて書いていますが、拡張されたユークリッドの互除法を文字式でやっているだけです。
そんでもって、bx＋ry＝1より、yについて解けばy＝(1－bx)／rになりますから、Y(b, r)＝(1－bX(b, r))／rが成り立つとのことなので、ax＋by＝1の解xは、x＝Y(b, r)＝(1－bX(b, r))／rと書けるというわけなのです。
xはX(a, b)と書けますから、X(a, b)＝(1－bX(b, r))／rという再帰関数になるのですね。rはa÷bのあまりでしたから、r＝a mod bと表わして、X(a, b)をf(a, b)と改めれば
f(a, b)＝(1－bf(b, a mod b))／(a mod b)が導けたことになります。
これだけだと再帰が止まりませんから、終端の値も求めなければなりません。
再帰していくといずれb＝1になるわけですが、b＝1になったら不定方程式はax＋y＝1となります。aが何であろうと(x, y)＝(0, 1)にしてしまえば特殊解になりますよね。ですから、b＝1のときはf(a, b)＝0と定義してやればいいことになります。
従って、冒頭の再帰関数の式と、特殊解の式が導かれます。yはax＋by＝1をyについて解いたものです。
ax＋by＝1の特殊解が分かっていれば、ax＋by＝cの特殊解は簡単に求まります。c倍すればいいだけですね。従って、ax＋by＝1の特殊解が(X, Y)であれば、ax＋by＝cの特殊解は(cX, cY)になります。
特殊解が分かれば、一般整数解もただちに求まります。
整数mを用いて、x＝cX－bm, y＝cY＋amと表わすことができますね。
最後に、一般整数解を用いて、一般自然数解を導きます。x, yが自然数であることからx＞0, y＞0が成り立ちますので、cX－bm＞0, cY＋am＞0が成り立つことになります。これをmについて解けば－cY／a＜m＜cX／bが導かれます。
床関数をf[ ]、天井関数をc[ ]と表わすと、その定義から
－cY／a＜f[－cY／a]＋1≦－cY／a＋1
cX／b－1≦c[cX／b]－1＜cX／b
がそれぞれ成り立ちますので、mの最小値、最大値は、それぞれ
f[－cY／a]＋1, c[cX／b]－1となりますから、
f[－cY／a]＋1≦m≦c[cX／b]－1としてmを整数の閉区間で表わすことができます。
このままでもいいのですけれど、なんとなく入れる整数は0, 1, 2, ……みたいにしたいので、始端をシフトして0から始まるように変数変換してやります。
全辺f[－cY／a]＋1を引いてやれば、
0≦m－f[－cY／a]－1≦c[cX／b]－f[－cY／a]－2
となりまして、間の変数をnと改めれば（n＝m－f[－cY／a]－1とすれば）、
0≦n≦c[cX／b]－f[－cY／a]－2という範囲ができあがります。
これならnの個数はc[cX／b]－f[－cY／a]－1で得られることは一目瞭然ですね＾＾；
これをさっき求めた一般整数解に代入してやることにより、次の一般自然数解を得られます。
これで解の公式が完成しました。
