RPGツクールと数学のブログ

RPGツクールと数学についてのブログです。

ベズー方程式の解の公式 その2(多変数 自然数解のみ)

自然数a, b, cからなる不定方程式ax+by=cの自然数解(x, y)を求める解の公式を前回作りました。

fermiumbay13.hatenablog.com

ここではこれを拡張して、任意の個数の元での自然数解を求めます。

ã¤ã¡ã¼ã¸ 1

k個の正の整数a[1], a[2], …, a[k](これらはそれぞれ互いに素であるとします)と、k個の変数x[1], x[2], …, x[k]との積和が、ある正の整数b[1]になるとき、このk個の変数がすべて自然数となるようなペアを求めることを考えます。

k=2だと普通のax+by=cの形の方程式ですね。k≧2とします。右辺の定数はb[1]としていますが、これは後の利便性のためです。

考え方としては次のように、それぞれを別の変数に置き換えることをして、ax+by=cの形に帰着させることにより解くということとします。

a[1]x[1]+a[2]x[2]+a[3]x[3]+a[4]x[4]+a[5]x[5]=b[1]

後ろの2項をb[4]と置く

a[4]x[4]+a[5]x[5]=b[4]
a[1]x[1]+a[2]x[2]+a[3]x[3]+b[4]=b[1]

さらに後ろの2項をb[3]と置く

a[4]x[4]+a[5]x[5]=b[4]
a[3]x[3]+b[4]=b[3]
a[1]x[1]+a[2]x[2]+b[3]=b[1]

さらに後ろの2項をb[2]と置く

a[4]x[4]+a[5]x[5]=b[4] … ④
a[3]x[3]+b[4]=b[3] … ③
a[2]x[2]+b[3]=b[2] … ②
a[1]x[1]+b[2]=b[1] … ①

これで、4本の式ができました。

順番が逆になってしまいましたが、右辺の定数をb[1]とおくと、上のようにきれいに式の添え字がそろうという仕組みなのですね。それぞれの方程式は前回の解の公式を用いて一般解を求めることが可能です。
①からx[1]とb[2]を求め、b[2]を用いて
②からx[2]とb[3]を求め、b[3]を用いて
③からx[3]とb[4]を求め、b[4]を用いて
④からx[4]とx[5]を求めて、すべての解を得る、というステップを踏みます。

このとき、それぞれ複数解がありますから、求める過程は入れ子構造となります。

①で求まった複数のb[2]に対してすべて②を適用、
②でも複数b[3]が求まるから、それらすべてに対して③を適用、……
という調子です。イメージできますよね。

ということで、このような働きをする解の公式を作りました。構造が上のようになっているだけで、実際に適用するときは代入するだけでOKです。

ちょっと前回と記号を変えています。

解の公式

ax+by=1の解を与える再帰関数
(画像ミスってて申し訳ないのですが、X0の右辺にあるfはX0のことです)
ã¤ã¡ã¼ã¸ 2

任意の整数n[1], n[2], …, n[k-1]
ã¤ã¡ã¼ã¸ 3

新たにおく定数b[i]の漸化式(選択する整数により決定)
ã¤ã¡ã¼ã¸ 4

求める一般自然数
ã¤ã¡ã¼ã¸ 5

こんなの手計算でやってられませんね。いちおう例をやってみましょう。

例: 23x+39y+55z+73w=500

上記の解の公式は、大きい係数のものを前にしたほうが計算しやすい構造になってます。そのため、ひっくり返したものとして計算しましょう。

それぞれの文字を次のように置いて、解の公式に代入します。

ã¤ã¡ã¼ã¸ 6

それでは順番に代入してみます。まずは最初の整数n[1]の範囲を求めます。

ã¤ã¡ã¼ã¸ 7

よって、0≦n[1]≦5となります。ひいい、6個もあるのか……

それぞれの場合において、b[2]を求めます。めんどくさいから同じ文字でb[2]としていますけれど、本当は選ぶ整数によって変わるので、b[2](0), b[2](1), … b[2](5)のように、関数みたいに表わした方がいいでしょうね。b[3]は2変数。b[3](2, 3)みたいな感じ。なんでもありません。

b[2]=62+73n[1]となるので、n[1]=0~5を順に代入してb[2]をそれぞれ求めます。

ã¤ã¡ã¼ã¸ 8

6つあるb[2]をそれぞれ用いて、整数n[2]の範囲をそれぞれ求めます。このとき、範囲の上限値が-1になったら解なしであることを意味します。これで枝切りします。

ã¤ã¡ã¼ã¸ 9

n[2]の範囲は、b[2]によって上のように分岐します。ヤバイですね。思ってたより分岐が多いです。なんで変数4個にしちゃったんだろう。大丈夫かなぁ……

続いてb[3]です。これはb[2]とn[2]によって分岐します。

ã¤ã¡ã¼ã¸ 10

n[2]を範囲内のぶんだけ代入していけばb[3]が求まりますね……

b[2]=62のとき
┗━━ n[2]=0のとき
    ┗━━ b[3]=7
b[2]=135のとき
┣━━ n[2]=0のとき
┃   ┗━━ b[3]=25
┗━━ n[2]=1のとき
    ┗━━ b[3]=80
b[2]=208のとき
┣━━ n[2]=0のとき
┃   ┗━━ b[3]=43
┣━━ n[2]=1のとき
┃   ┗━━ b[3]=98
┗━━ n[2]=2のとき
    ┗━━ b[3]=153
b[2]=281のとき
┣━━ n[2]=0のとき
┃   ┗━━ b[3]=6
┣━━ n[2]=1のとき
┃   ┗━━ b[3]=61
┣━━ n[2]=2のとき
┃   ┗━━ b[3]=116
┣━━ n[2]=3のとき
┃   ┗━━ b[3]=171
┗━━ n[2]=4のとき
    ┗━━ b[3]=226
b[2]=354のとき
┣━━ n[2]=0のとき
┃   ┗━━ b[3]=24
┣━━ n[2]=1のとき
┃   ┗━━ b[3]=79
┣━━ n[2]=2のとき
┃   ┗━━ b[3]=134
┣━━ n[2]=3のとき
┃   ┗━━ b[3]=189
┣━━ n[2]=4のとき
┃   ┗━━ b[3]=244
┗━━ n[2]=5のとき
    ┗━━ b[3]=299
b[2]=427のとき
┣━━ n[2]=0のとき
┃   ┗━━ b[3]=42
┣━━ n[2]=1のとき
┃   ┗━━ b[3]=97
┣━━ n[2]=2のとき
┃   ┗━━ b[3]=152
┣━━ n[2]=3のとき
┃   ┗━━ b[3]=207
┣━━ n[2]=4のとき
┃   ┗━━ b[3]=262
┣━━ n[2]=5のとき
┃   ┗━━ b[3]=317
┗━━ n[2]=6のとき
    ┗━━ b[3]=372

今度から変数は3個にしますね……

今回はk=4ですから、i=3となるときの整数n[3]の範囲は公式にあるように、分岐してより複雑になります。そこにそれぞれ代入してn[3]の範囲を求めるのですが、ここはもう省略してもいいですよね。

結論から言えば、b[2]=427のときのn[2]=4, 5, 6のとき以外はすべて解なしになります。
n[2]=4, 5, 6のとき、すべてn[3]=0となって、それぞれに解が1個ずつあることがわかります。
ここまでですべての整数の範囲が確定しますので、解の個数は3個であることが分かります。

後はもう道具は出そろっているので、解x[i]を求めていけばよいですね。
b[1]=500, n[1]=5, b[2]=427までは共通しているので、x[1]は全部同じです。
公式にあるように、x[1], x[2]は一番シンプルな式を使って求めます。

ã¤ã¡ã¼ã¸ 11

x[2]を求めるにはn[2]が必要なので、n[2]=4, 5, 6をそれぞれ代入しましょう。

ã¤ã¡ã¼ã¸ 12

最後の2個であるx[3], x[4]はめんどくさい式ですね。ここで初めてユークリッドの互除法を使ってます。ユークリッドの互除法を行う再帰関数の計算はa[3], a[4]だけで行えるので共通です。

ã¤ã¡ã¼ã¸ 13

X0が-10、Y0が17になりました。

ã¤ã¡ã¼ã¸ 14

あとはb[3], n[3]のペアを代入して、具体的な解を得ます。

ã¤ã¡ã¼ã¸ 15

ã¤ã¡ã¼ã¸ 16

長い道のりでしたね。これですべての解x[1]~x[4]が得られました。

何やってたのか忘れちゃってましたけど、23x+39y+55z+73w=500を解いていたのでしたね。もともとのはx[4], x[3], x[2], x[1]がx, y, z, wになるようにひっくり返していたので、解は次のようになります。

ã¤ã¡ã¼ã¸ 17

シンプルな解の割に合わない計算量です。

元が1個増えると再帰回数が1個増えるので、手計算ではこのぐらいが限界でしょう。コンピュータでプログラミングするには実装しやすいんじゃないですかね。

導出

Σ[i=1→k]a[i]x[i]=b[1](k≧3)を
a[i]x[i]+b[i+1]=b[i](1≦i≦k-2) … ①
a[k-1]x[k-1]+a[k]x[k]=b[k-1](i=k-1) … ②
というk-1本の式に分けます。冒頭に述べた通りですね。

ここで、ax+by=cの一般解を範囲が0≦n≦N(a, b, c)であるような整数nを用いて
x=X(a, b, c, n), y=Y(a, b, c, n)と表現します。

前回示した通り、a, bが互いに素でa, b, cが正の整数であるとき、
N(a, b, c)=c[cX0(a, b)/b]-f[-cY0(a, b)/a]-2
X(a, b, c, n)=cX0(a, b)-b(n+f[-cY0(a, b)/a]+1)
Y(a, b, c, n)=cY0(a, b)+a(n+f[-cY0(a, b)/a]+1)
X0(a, b)=(1-bX0(b, a mod b))/(a mod b)(b≠1), 0(b=1)
Y0(a, b)=(1-aX0(a, b))/b
としてそれぞれ表わすことができるのでした。f[ ]は床関数、c[ ]は天井関数です。ややこしいですね。

①, ②の一般解をそれぞれこの表現を用いて表わしてみましょう。

① a[i]x[i]+b[i+1]=b[i]

これは整数n[i]を用いて表現します。n[i]の範囲は0≦n[i]≦N(a[i], 1, b[i])と表わせて、一般解はx=X(a[i], 1, b[i], n[i]), y=Y(a[i], 1, b[i], n[i])と表わせます。

b[i+1]の係数が1であることを利用すると、これらは割と簡潔に表わせます。
X0(a[i], 1)=0
Y0(a[i], 1)=1
N(a[i], 1, b[i])=-f[-b[i]/a[i]]-2=c[b[i]/a[i]]-2
X(a[i], 1, b[i], n[i])=-b(n[i]+f[-b[i]/a[i]]+1)
Y(a[i], 1, b[i], n[i])=b[i]+a[i](n[i]+f[-b[i]/a[i]]+1)
ユークリッドの互除法を必要とする部分が必要なくなるので、式に代入するだけでそれぞれの解を得られるようになります。従いまして、
0≦n[i]≦c[b[i]/a[i]]-2
x[i]=-b(n[i]+f[-b[i]/a[i]]+1)
b[i+1]=b[i]-a[i](c[b[i]/a[i]]-1-n[i])
と求まります。b[i+1]の所はa[i]の係数をこっそり逆にしてます。これらが1≦i≦k-2のときの各値ですね。

② a[k-1]x[k-1]+a[k]x[k]=b[k-1]

整数n[k-1]の範囲は0≦n[k-1]≦N(a[k-1], a[k], b[k-1])となり、一般解は、x=X(a[k-1], a[k], b[k-1], n[k-1]), y=Y(a[k-1], a[k], b[k-1], n[k-1])と表わせます。

これらは単純に代入しただけで、そのおかげでずいぶん汚い式になっています。

0≦n[k-1]≦c[b[k-1]X0(a[k-1], a[k])/a[k]]-f[-b[k-1]Y0(a[k-1], a[k])/a[k-1]]-2
x[k-1]=b[k-1]X0(a[k-1], a[k])-a[k](n[k-1]+f[-b[k-1]Y0(a[k-1], a[k])/a[k-1]]+1)
x[k]=b[k-1]Y0(a[k-1], a[k])+a[k-1](n[k-1]+f[-b[k-1]Y0(a[k-1], a[k])/a[k-1]]+1)
これらがi=k-1のときの式です。

以上ですべて求められますね。

一般的な部分である1≦i≦k-2の方が、特別な部分であるi=k-1よりもシンプルな式になるのは違和感がありますね。