Shironetsu Blog

@shironetsuのブログ

フルーツ・パズルをSageMathで解く――楕円曲線論初歩

これは何

「95%の人が解けない問題」とされるインターネットミームに書かれた方程式は、見た目に反して背景に楕円曲線がある難問で、この記事ではSageMathという数式処理システムを利用して解いてみた。

SageMath Notebookを以下のページに上げています。

github.com


「95%の人が解けない問題」


https://qph.fs.quoracdn.net/main-qimg-5b0690e302a38cf2a8068158199e7a21-c
Figure 1. 95% of people cannot solve this! <equation> Can you find positive whole values for :apple:, :banana:, and :pineapple:?

記事タイトルのフルーツ・パズルとはFig. 1で問われているパズルのこと。数年前にミーム化した画像なので知っている人も多いだろう。自分はTwitterで知ったはず。要は
\begin{align}
\frac{X}{Y+Z} + \frac{Y}{X+Z} + \frac{Z}{X+Y} = 4
\end{align}の正整数解 {(X,\,Y,\,Z)} を見つけよ、という問題で、煽り気味に「95%の人が解けない!」という文句が付いている。しかし、実のところ見た目に反して「頭の体操」で済む難度の方程式ではない。実際、最小の解は
\begin{align}
X &= 4373612677928697257861252602371390152816537558161613618621437993378423467772036 \simeq 4.4 \times 10^{78},\\
Y &= 36875131794129999827197811565225474825492979968971970996283137471637224634055579 \simeq 3.7 \times 10^{79},\\
Z &= 154476802108746166441951315019919837485664325669565431700026634898253202035277999 \simeq 1.5 \times 10^{80}
\end{align}
と非常に巨大な数になる。

こんな大きな整数が一体どうやって見つかるのだろうか? {X,\,Y,\,Z} の3重ループで探索すれば、整数解としては {(X,Y,Z)=(-1,4,11), (-5,9,11)} などが見つかる。しかし正の解にたどり着くには、その前に太陽の寿命が尽きてしまう。

このフルーツ・パズルには既に優れた解説がいくつもある[1-5]。特に論文”"An unusual cubic representation problem"[3]には、より一般的な形の {a/(b+c)+b/(c+a)+c/(a+b)=N} の整数解に関する包括的な議論がある。

鍵は楕円曲線だ。いずれの解説でも、楕円曲線上の有理点を計算する問題へと翻訳するという点で、求解の手順は基本的に一致している。概要を述べよう。


(1)
\begin{align}
f(X,Y,Z) = \frac{X}{Y+Z} + \frac{Y}{X+Z} + \frac{Z}{X+Y} - 4
\end{align}とする。探しているのは {f(X,Y,Z)=0} の正整数解だが、任意の {t\neq 0}
\begin{align}
f(tX,tY,tZ) = f(X,Y,Z)
\end{align}だから、{(X,Y,Z)} が全て同符号の有理数解があれば適当な {t} (分子の最小公倍数)倍して整数解を得ることができる。

(2)
\begin{align}
g(X,Y,Z) = (X+Y)(Y+Z)(Z+X)f(X,Y,Z)
\end{align}として分母を払うと、{C:\,g(X,Y,Z)=0}射影平面 {\mathbb{P}^2} 内の三次曲線である。双有理変換 {\phi: (X,Y,Z)\mapsto(X',Y',Z')} を適切に選ぶと
\begin{align}
h(X',Y',Z') = h(\phi(X,Y,Z)) = g(X,Y,Z)
\end{align}で定義される {E:\,h(X',Y',Z')=0}楕円曲線の標準形になる。{E} の 有理点 {(X',Y',Z')} の逆像 {\phi^{-1} (X',Y',Z')} から {C} の 有理点が得られるから、前者を求める。

(3){E} の有理点(特に {E} 上の Mordel-Weil 群の生成元)を一つ見つけて {P}とする。楕円曲線上の加法によって 次々に{n} 倍点 {n P} を計算し、{C} の有理点に戻したとき全て同符号になるまで続ける。ネタバレすると {n=9} で上の解が見つかる。

――

以上が解を求めるための手順となる。種が分かっていれば方針は簡明とも言える。ただ、アルゴリズムが知られているとはいえ

  • 双有理変換 {\phi} の決定
  • 楕円曲線上の点の和の計算

は骨の折れる作業だ。kaityo256氏によるQiitaの記事 [5]では、双有理変換はMathematicaと手で、加法はフルスクラッチRubyコードの実装を行っている実例を見ることができる。

もっと簡単にこの求解の手順を体験できないだろうか? SageMath を利用すればそれができる。

www.sagemath.org

SageMath は オープンソース・無償の数式処理システム で、Python3ベースの文法が採用されている。

安定版の9.3(2021年5月25日現在)を上記公式サイトからダウンロードできる。また、クラウドサービスCoCalc ではブラウザ上、オンラインでSageMathの機能が利用できる。

cocalc.com

おすすめはJupyter Notebook上でSageMathが扱えるSageMath Notebookだ。この記事の内容は以下のファイルで試すことができる。

github.com

ではやってみよう。


SageMathで楕円曲線

まず {f(X,Y,Z)} を定義する。

入力:

f(X,Y,Z) = X/(Y+Z) + Y/(Z+X) + Z/(X+Y) - 4

ちなみにshow()を使うとlatex形式で式を表示できる。

入力:

show(f(X,Y,Z))

出力:
{
\begin{align}
\frac{X}{Y + Z} + \frac{Y}{X + Z} + \frac{Z}{X + Y} - 4
\end{align}
}

試しに単純な3重ループで整数解を探してみる。

入力:

N = 20
for X in range(-N,N+1):
    for Y in range(X,N+1):
        for Z in range(Y,N+1):
            if(X+Y==0 or Y+Z==0 or Z+X==0):
                continue
            if(f(X,Y,Z)==0):
                print((X,Y,Z))

出力:

(-11, -9, 5)
(-11, -4, 1)
(-5, 9, 11)
(-1, 4, 11)

確かにこれらは解になっている。これだけ見るとずいぶん穏当な問題に見えるが罠だ。

次に、{g(X,Y,Z)} を定義する。

入力:

g(X,Y,Z) = f(X,Y,Z)*(X+Y)*(Y+Z)*(Z+X)
g(X,Y,Z) = g(X,Y,Z).simplify_full().expand() #整理して展開
show(g(X,Y,Z))

出力:
{
\begin{align}
X^{3} - 3 \, X^{2} Y - 3 \, X Y^{2} + Y^{3} - 3 \, X^{2} Z - 5 \, X Y Z - 3 \, Y^{2} Z - 3 \, X Z^{2} - 3 \, Y Z^{2} + Z^{3}
\end{align}
}

探しているのはこれによって定義される {\mathbb{Q}} 上の3次曲線 {C:\,g(X,Y,Z)=0} の正の整数点だ。そのためには、{g(X,Y,Z)}{\mathbb{Q}} 係数の多項式であることを宣言する必要がある。

{\mathbb{Q}} 上の3変数 {X,Y,Z} による多項式環 R を以下のように定義する。

入力:

R = PolynomialRing(QQ, 'X,Y,Z'); R

出力:

Multivariate Polynomial Ring in X, Y, Z over Rational Field

QQ はSageMathにおける有理数体を表す予約語となっている。

R をコンストラクタとして先の g を渡すと3変数多項式 R が定義される。

入力:

cubic = R(g)

出力:

X^3 - 3*X^2*Y - 3*X*Y^2 + Y^3 - 3*X^2*Z - 5*X*Y*Z - 3*Y^2*Z - 3*X*Z^2 - 3*Y*Z^2 + Z^3

式としては同じものだが、gcubic とでは属す集合(と型)が異なる。parent() によってこれを見ることができる:

入力:

print(g.parent())
print(cubic.parent())

出力:

g:  Callable function ring with arguments (X, Y, Z)
cubic:  Multivariate Polynomial Ring in X, Y, Z over Rational Field

g は「ただの関数」だが、cubic は明示的に有理数体上の多項式として定義されている。

cubicから楕円曲線を作ろう。SageMathで楕円曲線を初期化する方法は色々ある。
doc.sagemath.org

例を挙げると、

  • EllipticCurve([a1,a2,a3,a4,a6])

「長い」Weierstrass標準形
\begin{align}
y^2+a_1xy+a_3y = x^3+a_2x^2 + a_4 x + a_6
\end{align}の係数によって定義。

  • EllipticCurve_from_j(j0)

{j}-不変量の値から定義。

  • EllipticCurve('label')

Cremona databaseと呼ばれる楕円曲線のデータベース上に登録された「ラベル」("11a"や"37b2"のような)から定義。

ちなみに、本記事で扱う楕円曲線Eは"910c1"である。E.cremona_label() から得られる。LMFDBというウェブサイトに詳細な情報が集約されていて、対応するページ(910c1はここ[8])を開くためのメソッドも存在する:E.lmfdb_page()

など。

今の場合、3次曲線から構成するので、読んで字の如くのEllipticCurve_from_cubicを利用する。

デフォルトではこの関数からは戻り値として双有理写像が与えられる。上で定義した {\phi} である。

入力:

phi = EllipticCurve_from_cubic(cubic); phi

出力:

Scheme morphism:
  From: Projective Plane Curve over Rational Field defined by X^3 - 3*X^2*Y - 3*X*Y^2 + Y^3 - 3*X^2*Z - 5*X*Y*Z - 3*Y^2*Z - 3*X*Z^2 - 3*Y*Z^2 + Z^3
  To:   Elliptic Curve defined by y^2 + x*y = x^3 + 69*x^2 + 1365*x + 8281 over Rational Field
  Defn: Defined on coordinates by sending (X : Y : Z) to
        (-X - Z : X : 6/91*X - 1/91*Y + 6/91*Z)

{C: g(X,Y,Z)=0} が非同次形({\phi} の像 {(X',\,Y',\,Z')} に対して {x=X'/Z',\, y=Y'/Z'} と置き換え)で
\begin{align}
E:\,y^2 + xy = x^3 + 69x^2 + 1365x + 8281
\end{align}という楕円曲線と双有理同値であると言っている。また、{\phi} が具体的に
\begin{align}
\phi:\,\left\lbrack X:\, Y:\, Z\right\rbrack \mapsto \left\lbrack -X-Z:\, X:\, \frac{6}{91}X-\frac{1}{91}Y+\frac{6}{91}Z\right\rbrack
\end{align}という形で与えられている。つまり、
\begin{align}
x &= \frac{91(-X-Z)}{6X-Y+6Z},\\
y &= \frac{91X}{6X-Y+6Z}
\end{align}
とすれば楕円曲線 {E} は元の {C} に一致するはずである。詳細は省くが、以下の操作でこれを確かめられる:

入力:

E.defining_polynomial()(phi.defining_polynomials()) * phi.post_rescaling()


出力:

X^3 - 3*X^2*Y - 3*X*Y^2 + Y^3 - 3*X^2*Z - 5*X*Y*Z - 3*Y^2*Z - 3*X*Z^2 - 3*Y*Z^2 + Z^3

{E} の「定義多項式」(上で定義した {h})に {\phi} の像を代入した結果が元の {g(X,Y,Z)} に帰ってきている。

さて、{E} の有理点を探そう。Mordell-Weil 群の生成元である有理点が gens() から得られる。

入力:

E.gens()

出力:

[(-39 : 52 : 1)]

今の場合位数無限の生成元は一つしかない。このことは {E}ランク{1} であることを反映している:

入力:

E.rank()

出力:

1

そこで、この点を改めて {P} とおく。

入力:

[ P ] = E.gens()

{P}{\phi} の逆写像 {\phi^{-1}} によって引き戻せば {C} 上の有理点になるはずだ。これは実際、以下のように計算できる:

入力:

phi_inv = phi.inverse() #逆写像
phi_inv(P)

出力:

(-4 : -11 : 1)

{(X,\,Y,\,Z)=(-4,-11,1)}{C} 上の有理点になることは3重ループによる計算で見た通りだった。

楕円曲線上の加法によって {P+P=2P} で与えられる点とその逆像は、

入力:

print(2*P)
print(phi_inv(2*P))

出力:

(1859/25 : -123487/125 : 1)
(-9499/8784 : -5165/8784 : 1)

と計算される。後者を {8784} 倍して {(-9499,\, -5165,\, 8784)} とするとこれは {f(X,Y,Z)=0} の整数解になっている。この段階ですでに素朴なループではhoursのオーダーを要する領域に入っている。

同様に、{P}{n} 倍点 {P+P+\cdots+P=nP}の逆像から {C} の有理点が次々に得られる。ループで回そう。

入力:

for n in range(10):
    a, b, c = phi_inv(n*P)
    k = lcm(a.denominator(), b.denominator())
    p, q, r = k*a, k*b, k*c  # c=1 で規格化されているため、aとbの最小公倍数kをかける
    if((sgn(p), sgn(q), sgn(r)) == (-1, -1, -1)): # 同符号だが負数のとき(-1)倍する
        p, q, r= -p, -q, -r
    if((sgn(p), sgn(q), sgn(r)) == (1, 1, 1)):
        print('>>>positive solution !<<<')
    print(n, (p, q, r))

出力:

0 (-1, 0, 1)
1 (-4, -11, 1)
2 (-9499, -5165, 8784)
3 (-679733219, -883659076, 375326521)
4 (-6696085890501216, -6334630576523495, 6531563383962071)
5 (-5824662475191962424632819, -5048384306267455380784631, 2798662276711559924688956)
6 (-287663048897224554337446918344405429, -434021404091091140782000234591618320, 399866258624438737232493646244383709)
7 (-3386928246329327259763849184510185031406211324804, -1637627722378544613543242758851617912968156867151, 678266970930133923578916161648350398206354101381)
8 (-343258303254635343211175484588572430575289938927656972201563791, -2110760649231325855047088974560468667532616164397520142622104465, 2054217703980198940765993621567260834791816664149006217306067776)
>>>positive solution !<<<
9 (154476802108746166441951315019919837485664325669565431700026634898253202035277999, 4373612677928697257861252602371390152816537558161613618621437993378423467772036, 36875131794129999827197811565225474825492979968971970996283137471637224634055579)

こうして晴れて {n=9} で正整数解が見つかった。ちなみに次の解は {n=17} で、

(1440354387400113353318275132419054375891245413681864837390427511212805748408072838847944629793120889446685643108530381465382074956451566809039119353657601240377236701038904980199109550001860607309184336719930229935342817546146083848277758428344831968440238907935894338978800768226766379, 
1054210182683112310528012408530531909717229064191793536540847847817849001214642792626066010344383473173101972948978951703027097154519698536728956323881063669558925110120619283730835864056709609662983759100063333396875182094245046315497525532634764115913236450532733839386139526489824351, 
9391500403903773267688655787670246245493629218171544262747638036518222364768797479813561509116827252710188014736501391120827705790025300419608858224262849244058466770043809014864245428958116544162335497194996709759345801074510016208346248254582570123358164225821298549533282498545808644)

である。{XYZ \simeq 10^{856}}


ねじれ部分

ところで {X,Y,Z} の入れ替えも明らかに {f(X,Y,Z)=0} の整数解になる。これはどこに現れるのだろう?

{E} のMordell-Weil 群の生成元は {P} だけではない。{E} にはねじれ部分もあって、以下のように得られる。

入力:

E.torsion_points()

出力:

[(-14 : 7 : 1),
 (-13 : 0 : 1),
 (-13 : 13 : 1),
 (0 : -91 : 1),
 (0 : 1 : 0),
 (0 : 91 : 1)]

6個のねじれ点が得られた。すなわち {E} の Mordell-Weil群は {\mathbb{Z}\times \mathbb{Z}/6\mathbb{Z}} と同型ということになる。なお、この6という値は E.torsion_order() からも出せる。

ねじれ点の位数(すべて有限である)はそれぞれ以下のように計算できて、生成元として位数 {6}(0 : 91 : 1)を取れる。

入力:

for T in E.torsion_points():
    print(T.order(), T)

出力:

2 (-14 : 7 : 1)
3 (-13 : 0 : 1)
3 (-13 : 13 : 1)
6 (0 : -91 : 1)
1 (0 : 1 : 0)
6 (0 : 91 : 1)

入力:

Q = E(0, 91, 1) #ねじれ部分群の生成元

{E} の Mordell-Weil群はこの2点 {P,\,Q} によって生成される。{\phi^{-1}(\pm P + n Q) } の逆像を見てみる。

入力:

for m in [-1, 1]:
    for n in range(E.torsion_order()):
        a, b, c = phi_inv(m * P + n * Q)
        k = lcm(a.denominator(), b.denominator())
        a, b, c = k*a, k*b, k*c
        print((m, n), (a,b,c))

出力:

(-1, 0) (-1, 11, 4)
(-1, 1) (-5, 9, 11)
(-1, 2) (4, -1, 11)
(-1, 3) (11, -5, 9)
(-1, 4) (-11, -4, 1)
(-1, 5) (-9, -11, 5)
(1, 0) (-4, -11, 1)
(1, 1) (-5, 11, 9)
(1, 2) (-1, 4, 11)
(1, 3) (9, -5, 11)
(1, 4) (11, -1, 4)
(1, 5) (-11, -9, 5)

\begin{align}
(m,n)=(-1,0),\, (-1,2)\, (-1,4),\, (1,0)\, (1,2),\, (1,4)
\end{align}が {(X,\,Y,\,Z)=(-1,11,4)} への置換(全体にかかる {(\pm 1)} は無視していい)に、
\begin{align}
(m,n)=(-1,1),\, (-1, 3)\, (-1,5),\, (1,1)\, (1,3),\, (1,5)
\end{align}が {(X,\,Y,\,Z)=(-5, 9, 11)} への置換に対応している。このことは、{E} 上の点 {S} への作用として、
\begin{align}
S &\mapsto S,\\
S &\mapsto S+2Q,\\
S & \mapsto S+4Q,\\
S &\mapsto -S,\\
S &\mapsto -S+2Q,\\
S &\mapsto -S+4Q
\end{align}
の6つの写像が合成で閉じていて({6Q = O} に注意)、{3} 次対称群と同型であることの表れとなっている。

{E} 上の加法と引き戻しで {f(X,Y,Z)}=0 の解に対する置換が引き起こされていることを具体的に示したかったが力及ばなかった。楕円曲線上の加法をシンボリックに行う方法をご存知の方がいれば教えてください。

なお、{\phi^{-1}(13P+Q)} から別の正整数解が得られる。

(184386514670723295219914666691038096275031765336404340516686430257803895506237580602582859039981257570380161221662398153794290821569045182385603418867509209632768359835, 
16666476865438449865846131095313531540647604679654766832109616387367203990642764342248100534807579493874453954854925352739900051220936419971671875594417036870073291371, 
32343421153825592353880655285224263330451946573450847101645239147091638517651250940206853612606768544181415355352136077327300271806129063833025389772729796460799697289)

がそれで、{XYZ\simeq 10^{500}} と、{ \phi^{-1}(9P),\, \phi^{-1}(17P)} それぞれから得られる解の間の大きさになる。


フルーツ・パズルの発祥はどこ?

ところで、フルーツ・パズルはどこからやってきたのだろうか。

2019年3月のRedditのスレッド "I'm stuck on this one" [6]にてこのフルーツ・パズルが議論されていて、以下の投稿がある:

このTCGSonIce氏による投稿で発祥地とされているのが、2017年1月の以下のスレッド "[request/fun] I'm really sick of all the facebook fruit math bull that's going on lately. Does anyone want to create a truly difficult math problem with pictures of fruit to counter this?"[7] 。スレッドタイトルはだいたい「最近フェイスブックで流行っている、くだらないフルーツ数学クイズには心底うんざりしている。対抗して、フルーツの絵を使った本当に難しい数学の問題を誰か作らないか?」といったところ。要は釣りスレだ。どうもこのころのフェイスブックにそういう流れがあったらしい。

このスレッドに以下の投稿がある。

Obyeag氏が2016年1月のmathoverflowの質問と回答 [2]を元に問題を考案し、louiswins氏によってa,b,cをリンゴ・バナナ・パイナップルに置き換えた画像がアップロードされている。

View post on imgur.com
imgur.com
なるほど解像度が高い。この段階では正整数解という条件は画像には直接書き込まれておらず、説明文に付記したとされている。数式以外明らかにフォントが違うので、拡散される中でどこかのまた別の誰かが書き換えたのだろう。それがどこかなのかまでは分からなかった。ちなみに(実際ネイティブがどう受け取るのかは分からないが、)"positive whole values"というのはふつう学術的には使われない表現のはずで、うまい隠蔽になっているように思う。

つまりまあミーム兵器として製造されたのがこのフルーツ・パズルだったわけだ。我々とこの素敵な方程式との出会いの裏には、こんな奇妙な経緯があったのだった。


References