この記事を書いたきっかけ

サイモン・シン著の『数学者たちの楽園　「ザ・シンプソンズ」を作った天才たち』（新潮社）を読んで思い立ちました．
登場した数式のひとつが面白くて，自分でも見つけてみたいなあと思ったことがきっかけです．
この本の内容は”シンプソンズ”の脚本家チームには高度な数学のバックグラウンドを持つメンバーが多く在籍し，彼らは作品の中に数多くの数学ジョークを仕込んでいたというものです．めちゃくちゃ面白かったのでオススメです．

ホーマーの最終定理?

”シンプソンズ”の登場人物ホーマー・シンプソンはことあるごとに怪しげな発明をして周囲を巻き込みますが，<エバーグリーン・テラスの魔法使い>(1998)の回ではその発明熱が頂点に達します．
そして，この回でホーマーが黒板に書いた数式が衝撃的なものです．

特に注目したいのは上から2つめの数式．
$$ 3987^{12} + 4365^{12} = 4472^{12} $$
なんとこの式，フェルマーの最終定理を覆しています！
フェルマーの最終定理とは
$$ x^{n} + y^{n} = z^{n} \;(n>2) $$
を満たす整数の組$[x,y,z]$は存在しないというものです．

この定理の証明は1994年に数学者アンドリュー・ワイルズが発表していますが，ホーマーが書いた数式は一体何なのでしょうか？

ホーマーはフェルマーの最終定理の反例を見つけて，数学史に残る偉業を成し遂げてしまったのでしょうか？

（サイモン・シンにならって以下，本記事では上記の式を「フェルマー方程式」と呼びます）

電卓で確かめてみると？

等式$$ 3987^{12} + 4365^{12} = 4472^{12} $$が本当に正しいか，電卓で確かめてみます．

フェルマー方程式のニアミス解

そう，ここにカラクリが潜んでいます．
表示が12桁よりも多い電卓を使うと結果が変わってくるのです．次の式に現れる値に近くなります．
$$ 3987^{12} + 4365^{12} = 4472.0000000070576171875^{12} $$実は先のホーマーの式はフェルマーの方程式に対する「ニアミス解」なのです！

この組み合わせは脚本家のひとり，デーヴィッド・S・コーエンが整数を次々とチェックするコンピュータ・プログラムを書いて見つけたそうです．

なんとこの方，カリフォルニア大学バークレー校で大学院生をしていた頃に，ワイルズによるフェルマーの最終定理の証明で非常に重要な役割を果たしたケン・リベットの講義に出席していたとのこと．

羨ましい限りです！

Pythonでニアミス解を探す

さて，ここまで長くなりましたが，Python3で同様にフェルマー方程式のニアミス解となるような整数組を探すプログラムを書きます．
ちなみに『数学者たちの楽園』の巻末にはコーエンがニアミス解を探すためにC言語で書いたコードが掲載されています．
しかし，今回それはあまり読まず，お気に入りの言語を使って自分で書いてみます．

以下のような仕様にします．

• 整数nを入力として受け取る
• フェルマー方程式のニアミス解として小数点第n位の精度の整数組[x,y,z]を総当りで探す

ソースコード

以下のようになりました．愚直に総当りで求めます．
一応，次の２つの簡単な事実のみ使っています．

• $x < z,\; y < z$ であること
• $x,\;y$は対称性を持つため$x \leqq y$としても良いこと

今後，もっと計算を速くする工夫を入れようと思っています．

2018/03/18 更新

掲載していたコードを@shiracamus さんが綺麗に書き直してくれました．可読性・再利用性ともにかなり向上しています．感謝するとともに，ここに掲載させていただきます:)

#! -*- coding:utf-8 -*-
# This program finds a Fermat Number.
# You input an integer n at the first.

import sys
from itertools import count

def finder(n):
    e = 0.1**n
    return ((x, y, z, k, p)
            for z in count(3)
            for k in range(5, 12)
            for y in range(1, z)
            for x in range(1, y)
            for p in [(x**k + y**k) ** (1 / k)]
            if z < p
            if abs(z - p) < e)

def main():
    n = int(input('Input an integer n : '))
    x, y, z, k, p = next(finder(n))
    print('[x,y,z] =', [x, y, z])
    print('k =', k)
    print(x, '^', k, '+', y, '^', k, '=', p, '^', k, sep='')

if __name__ == '__main__':
    main()

実行結果

実際にニアミス解を求めてみます．

n = 5

まずは小数点第5位までの精度．

$ python3 Fermat_number.py 
Put an integer n : 5
[x,y,z] = [167, 192, 215]
k = 4
167^4+192^4=215.00000482976324^4

$$167^4+192^4=215.00000482976324^4$$
らしいです．
ニアミス解といえばニアミス解ですが，ちょっと粗いですね．

n = 6

続けて小数点第6位の精度で試します．

$ python3 Fermat_number.py
Put an integer n : 6
[x,y,z] = [419, 462, 477]
k = 10
419^10+462^10=477.00000087334587^10

$$419^{10}+462^{10}=477.00000087334587^{10}$$
という結果です．
10乗が登場して何となく嬉しいです．

n = 7

n=7で試します．かなり実行に時間がかかってしまいますが，なんとかニアミス解を得られます．

$ python3 Fermat_number.py
Put an integer n : 7
[x,y,z] = [1125, 2335, 2347]
k = 5
1125^5+2335^5=2347.0000000364907^5

$$1125^5+2335^5=2347.0000000364907^5$$
を得ました．

参考文献

サイモン・シン『数学者たちの楽園　「ザ・シンプソンズ」を作った天才たち』(青木薫 訳) 新潮社