はじめに
ちょっと前に、こんな問題が話題になった。
ネタ元はどこか知らないが、僕は以下のサイトで知った。
適当に和訳してみた。
元の出題は果物とか使ってて「数学パズルです〜」みたいな顔をしているが、これは典型的ディオファントス方程式で、途中で楕円曲線が出てくるガチな問題になっている。解説は先のサイトに書いてあるが、いくつか天下りな箇所があるので、そこをちょっと丁寧に解説しつつ、Rubyで解を探してみる。
方針
問題は以下の通り。
\frac{a}{b+c} +
\frac{b}{c+a} +
\frac{c}{a+b} = 4
上記を満たす「自然数」の組(a,b,c)を求める。すぐに分かる通り、もしひとつの解の組(a,b,c)が求まったら、それぞれ定数倍した(ta, tb, tc)も解となる。従って、大事なのは(a,b,c)の比であり、例えば$c=1$に固定して、(a,b)の有理数解を求める、という問題と等価である。従って、この問題の自由度(実質的な変数の数)は2となる。
さっきの式の分母を払って整理する。
a^3 + b^3 + c^3 -3(a^2b+ ab^2 +a^2c+ac^2 + b^2c+bc^2) - 5 abc =0
見て分かる通り、これは三次関数となる。二変数三次関数ということで、この時点で「楕円曲線が出てきそうだな」ということがわかる。
この式を適当に変形して、ワイエルシュトラスの標準形に持っていく。具体的には$(a,b,c)$を適当な線形変換で$(x,y,z)$に持っていき、$x/z$と$y/z$だけの式にすることで$z$を消去し、楕円曲線にする。
楕円曲線が得られたら、なんとかして一つ、有理数解を見つける。自然数解を見つけるのは極めて大変だが、有理数解まで広げれば一つ見つけるのはなんとかなる。もし、見つかった有理数解$(a,b,c)$が全て正であれば、分母を通分してやれば分子が自然数解を与える。しかし、すぐに見つかる解の組は正負混じった数になるため、そこから楕円曲線の性質を使って別の解を見つける。見つけた解が求める解でなければ、さらに別の解を・・・と続けて、最終的に求める解が見つかるまでこのプロセスを続ける、というのが解法の概要である。
標準形への変形
まず、ワイエルシュトラスの標準形へ変形する。ちゃんとした標準形にするのは面倒だが、ほどほどの形にするのはさほど難しくない。ポイントは、三変数を一気にやらず、二変数ずつ処理することである。
元々の式は３次元空間$(a,b,c)$の中の曲面を表していた。そこで、まず$c=0$を代入してしまう。これは曲面と$c=0$の平面と交差する曲線を表す。$(a,b)$の満たすべき式は、
\frac{a}{b} + \frac{b}{a} = 4
整理して、
a^2 - 4 a b + b^2 = 0
ここで、$(a,b)$を線形変換して項の数を減らす。対称的な形なので
\begin{align}
a &= x + y \\
b &= x - y
\end{align}
という変換を思いつくのは難しくないと思う。代入して整理すると、
x^2 = 3 y^2
項の数が3から2に減った。
さて、この変換を$c$を残した式に突っ込んで整理すると、
c^3 - 6 c^2 x - 11 c x^2 - 4 x^3 + (-c + 12 x) y^2 = 0
となる。ワイエルシュトラスの標準形の気持ちは、$x/z$と$y/z$だけの式に落とすことであったから、$x$と$y$が混ざる項が出てきては困る。うまい具合にそんな項は $(-c + 12 x) y^2$だけなので、
z = -c + 12 x
という変換を施せば、$x$と$y$が混ざる項が消える。面倒なのでMathematicaにやらせよう。
In[3]:= Expand[
c^3 - 6 c^2 x - 11 c x^2 - 4 x^3 + (-c + 12 x) y^2 /. c -> -z + 12 x]
Out[3]= 728 x^3 - 277 x^2 z + y^2 z + 30 x z^2 - z^3
すなわち、
$$
728 x^3 - 277 x^2 z + y^2 z + 30 x z^2 - z^3 = 0
$$
が得られた。これは$x,y,z$の同次形となっているから、両辺$z^3$で割って、$-x/z \equiv x$、$y/z \equiv y$として整理してやると1、
y^2 = 728 x^3 + 277 x^2 + 30 x + 1
なお、変換行列は以下の通り。
\begin{pmatrix}
a \\ b \\ c
\end{pmatrix}
=
\begin{pmatrix}
-1 & +1 & 0 \\
-1 & -1 & 0 \\
-12 & 0 & -1
\end{pmatrix}
\begin{pmatrix}
x \\ y \\ z
\end{pmatrix}
有理数解$(x,y,z)$が求まったら、この変換により$(a,b,c)$が求まる。有理数ベクトルの線形変換なので、やはり有理数ベクトルが得られる。もう少し頑張ればちゃんとしたワイエルシュトラスの標準形が得られるが、とりあえずこれで良いことにする。
有理数解探索
さて、得られた
y^2 = 728 x^3 + 277 x^2 + 30 x + 1
について、まずはなんでもいいから有理数解を求める必要がある。後の都合から、$x<0$の領域で探そう。まず探す領域を調べるため、$y=0$の場合の$x$の解を３つ調べる。ニュートン法でもgnuplotでも何使っても良いが、手元にMathematicaがあったのでそいつにやらせよう。
In[6]:= NSolve[728 x^3 + 277 x^2 + 30 x + 1 == 0, x]
Out[6]= {{x -> -0.22377}, {x -> -0.0798013}, {x -> -0.0769231}}
というわけで、$-0.22377 < x < -0.0798013$の領域で有理数解を探す。後の便利のために、以下のようなメソッドをNumericクラスとRationalクラスに追加しておこう。
class Numeric
def cuberoot
a = self
as = 1
ae = a
until as == ae or ae == as+1
an = (as+ae)/2
if an**3 < a
as = an
else
ae = an
end
end
return as if as**3 == a
return ae if ae**3 == a
return 0
end
def sqrt
a = self
as = 1
ae = a
until as == ae or ae == as+1
an = (as+ae)/2
if an**2 < a
as = an
else
ae = an
end
end
return as if as**2 == a
return ae if ae**2 == a
return 0
end
def is_square
return (self.sqrt**2 == self)
end
def is_cubic
return (self.sqrt**3 == self)
end
end
class Rational
def is_square
self.denominator.is_square && self.numerator.is_square
end
def is_cubic
self.denominator.is_sqrt && self.numerator.is_sqrt
end
def sqrt
d = self.denominator
n = self.numerator
Rational(n.sqrt,d.sqrt)
end
def cuberoot
d = self.denominator
n = self.numerator
Rational(n.cuberoot,d.cuberoot)
end
end
def abc(x,y,z)
a = -x +y
b = -x - y
c = -12 * x -z
t1 = a.denominator.lcm(b.denominator)
t = t1.lcm(c.denominator)
a = a*t
b = b*t
c = c*t
return a.numerator,b.numerator,c.numerator
end
超手抜きの有理数の平方根、立方根を返すメソッドである。ついでに、あとで使うので$(x,y,z)$を$(a,b,c)$に直すメソッドもここに追加しておく。
これを使って有理数解が探せる。先程の解が存在しそうなあたりを適当な有理数を突っ込んで探索するスクリプトはこんな感じ。
require './util.rb'
xmax = - 0.0769231
xmin = - 0.22377
def f(x)
728 * x**3 + 277 *x**2 + 30 *x + 1
end
(5..30).each do |d|
(1..(d-1)).each do |n|
f = - n.to_f/d.to_f
next if f < xmin
next if f > xmax
x = Rational(-n,d)
y = f(x)
next if !y.is_square
y = y.sqrt
a,b,c = abc(x,y,1)
puts "(#{a},#{b},#{c}) : #{x}, #{y}"
end
end
要するに$x$として適当な有理数を突っ込んでみて、それが有理数の平方数になっていれば$y$が求まるので、$(x,y,1)$を先程の行列に突っ込むことで$(a,b,c)$を求めるスクリプトになっている。ただし、最後に分母を払って整数解にしている。
実行するとこんな感じ。
$ ruby search.rb
(11,-5,9) : -1/9, 8/27
(9,-5,11) : -2/13, 7/13
(4,-1,11) : -3/14, 5/14
(11,-5,9) : -1/9, 8/27
(11,9,-5) : -2/25, 1/125
(9,-5,11) : -2/13, 7/13
(11,-5,9) : -1/9, 8/27
(4,-1,11) : -3/14, 5/14
ちゃんと求まっていますね。後で使うために$(x,y)$の値も表示している。
有理点生成
さて、ここからが問題である。楕円曲線には、「有理係数の楕円曲線上の有理数点を結ぶ直線と交わる３つめの交点も有理点である」という性質がある。解説は、例えばこのページを参照2。
先程、楕円曲線上でいくつか有理点が見つかったので、適当な二つを結んだ直線を引き、交点を既知の有理点リストに加え・・・ということを繰り返していき、どこかで対応する$(a,b,c)$が全て正(もしくは負)になる点が見つかれば終了、という手続きとなる。
とりあえず一つやってみよう。楕円曲線では、二つの有理点P1、P2を結ぶ直線と交わった点をP3'、そのy座標を反転させた点をP3とし、
$$
P1 + P2 = P3
$$
と表記する。また、P1とP2は同じ点であっても良い。この場合は点Pにおける接線と楕円曲線の交点を2P'とし、その座標を反転させた点を2Pとして
$$
P + P = 2P
$$
と表記する。表記上は加算となっているが、別に点のベクトルを足しているわけではないので注意。図示するとこんな感じ。
ここで、交点を$2P$とせず、x軸に対して反転した点を$2P$と呼ぶのは、次に使うため。
いま、$P$、$2P$という二つの有理点が求まったため、その二点を通る直線と楕円曲線の交点を$3P'$とし、その軸反転した点を$3P$とすると、また$P$と$3P$を使って別の有理点を作ることができる。
こうして次々に新しい有理点を構築して行くと、いつかは自然数解に対応する有理点に行き着くのでは？というのが元の問題の解き方となる。
さて、最初の点Pとして(-3/14, 5/14)を選ぼう。今、楕円曲線は
y^2 = 728 x^3 + 277 x^2 + 30 x + 1
なので、点(x,y)における接線の傾きは
\frac{dy}{dx} = \frac{1092 x^2 + 277 x + 15}{y}
である。点P(x_s, y_s)における楕円曲線の傾きを$l$とすると、点Pにおける接線は
y - y_s = l (x-x_s)
で表される。この式と楕円曲線から$y$を消去し、$x$の解を求めればよろしい。その為には三次方程式の解を求める必要がある。三次方程式の有理数解をちゃんと求めるのは面倒だが、ここでは点$P$で重解を持っていることがわかっているので、三次方程式は$(x+3/14)^2$の因子を持っていることがわかる。従って、これで多項式の因数分解をしてやれば、すぐに新たな解$2P'$が求まる。
こうして$2P'$が求まったら、軸反転して$2P$とし、次は$P$と$2P$を通る直線から$3P'$を求める。これも、既に$P$、$2P$のx座標という解がわかっているので、因数分解により三つ目の解がすぐに求まる。
図を見れば、新たに得た有理点を軸反転した点を次に使わないと、新しい解が見つからないのがわかると思う。$P$と$2P$を使って$3P'$を見つけた後、次に$P$と$3P'$を使ったとしたら、見つかるのが$2P$になるので。
スクリプトはこちら。
require './util.rb'
def diff(a3,a2,a1,a0,x,y)
Rational(3*a3*x**2 + 2 * a2*x + a1,2*y)
end
def calc(a3,a2,a1,a0,x)
y = a3*x**3+a2*x**2+a1*x+a0
y.sqrt
end
def findnext(a3,a2,a1,a0,x1,y1,x2,y2)
l = 0
if x1==x2
l = diff(a3,a2,a1,a0,x1,y1)
else
l = (y2-y1)/(x2-x1)
end
b3 = a3
b2 = a2 - l**2
b1 = a1 - 2*l*(y2-l*x2)
b0 = a0 - (y2-l*x2)**2
b2 = b2 + b3 * x1
b1 = b1 + b2 * x1
b0 = b0 + b1 * x1
b2 = b2 + b3 * x2
b1 = b1 + b2 * x2
nx = -b2/b3
ny = calc(a3,a2,a1,a0,nx)
return nx,ny
end
a3 = Rational(728,1)
a2 = Rational(277,1)
a1 = Rational(30,1)
a0 = Rational(1,1)
x = Rational(-3,14)
y = Rational(5,14)
nx = x
ny = y
9.times do |i|
a,b,c = abc(nx,ny,1)
puts "#{i+1}P: (#{a},#{b},#{c})"
nx, ny = findnext(a3,a2,a1,a0,x,y,nx,ny*(-1)**i)
end
実行するとこうなる(一瞬で終わります)。
1P: (4,-1,11)
2P: (8784,-9499,-5165)
3P: (679733219,-375326521,883659076)
4P: (6531563383962071,-6696085890501216,-6334630576523495)
5P: (5824662475191962424632819,-2798662276711559924688956,5048384306267455380784631)
6P: (399866258624438737232493646244383709,-287663048897224554337446918344405429,-434021404091091140782000234591618320)
7P: (3386928246329327259763849184510185031406211324804,-678266970930133923578916161648350398206354101381,1637627722378544613543242758851617912968156867151)
8P: (2054217703980198940765993621567260834791816664149006217306067776,-343258303254635343211175484588572430575289938927656972201563791,-2110760649231325855047088974560468667532616164397520142622104465)
9P: (154476802108746166441951315019919837485664325669565431700026634898253202035277999,36875131794129999827197811565225474825492979968971970996283137471637224634055579,4373612677928697257861252602371390152816537558161613618621437993378423467772036)
ちゃんと9Pで自然数の組を見つけることができた。ちなみに楕円曲線上の点の軌跡はこんな感じになる。
行ったり来たりしながら谷間に落ちていってる感じですね。
まとめ
ディオファントス方程式を、RubyやMathematicaを使って解いてみた。フェルマーの最終定理やコラッツ予想など、一見簡単そうに見える方程式が滅茶苦茶難しい、というのは整数論でよく見られる。こういうものをガチで解決するのは天才の仕事だが、そういうもので遊ぶのは我々にもできるし楽しい。整数論は紙と鉛筆だけでなく、計算機を使うともっと楽しめることが多いので、興味のある方はいろいろ試して見られたい。