ABC予想をPythonで調べる

もちろん宇宙際Teichmüller理論とか全然分からない。
計算機で何か分かればいいなというコンセプトでWikipediaをなぞりながら逐一プログラムを書いてゆく。

ABC予想 - Wikipedia
数学の難問「ABC予想」、京大教授が解明か  :日本経済新聞
ABC Conjecture « Programming Praxis

やってることは話しにならないお遊びなのでもっとかっちりとした説明を読みたい人は下リンクを参照。

根基とは何か

まずABC予想が何であるかを知る前に、根基という言葉の意味が何であるかを知る必要がある。
根基とは、

正の整数 n について、n の互いに異なる素因数の積を n の根基 (radical) と呼び、rad(n) と書く。

  • rad(16) = rad(2^4) = 2,
  • rad(17) = 17,
  • rad(18) = rad(2·3^2) = 2·3 = 6.
http://ja.wikipedia.org/wiki/ABC%E4%BA%88%E6%83%B3

である。
よって、根基を求めるためには、

  1. 素因数分解を行い、
  2. ユニークな素因数の積を求める

という手順を踏むことが分かる。
まず素因数分解を行うコードは以下のようになる。Python本家サイトのデモから拝借した。ナイーブな実装。

# -*- coding: utf-8 -*-
# http://svn.python.org/projects/python/trunk/Demo/scripts/fact.py
from math import sqrt
def factorize(n):
    if n < 1:
        raise ValueError('fact() argument should be >= 1')
    if n == 1:
        return []  # special case
    res = []
    # Treat even factors special, so we can use i += 2 later
    while n % 2 == 0:
        res.append(2)
        n //= 2
    # Try odd numbers up to sqrt(n)
    limit = sqrt(n+1)
    i = 3
    while i <= limit:
        if n % i == 0:
            res.append(i)
            n //= i
            limit = sqrt(n+1)
        else:
            i += 2
    if n != 1:
        res.append(n)
    return res

if __name__ == "__main__":
    print factorize(1000) # => [2, 2, 2, 5, 5, 5]

続いてrad(n)を実装する。これは素因数分解した数の配列から同じ素因数を除き積を求めればよい。

from operator import __mul__
def rad(n):
    return reduce(__mul__, set(factorize(n)))

これで根基を求めるコードは完成した。

ABC予想

続いてABC-tripleという定義が導入される。

自然数の三つ組 (a, b, c) で、a + b = c, a < bで、a と b は互いに素 (coprime) を満たすものを abc-triple と呼ぶ。

http://ja.wikipedia.org/wiki/ABC%E4%BA%88%E6%83%B3

互いに素、とはaとbが-1と1以外に約数を持たない場合のことを言う。互いに素であるかどうかを確かめるにはaとbの最大公約数を調べればよい。プログラミングにおいて最大公約数を求める方法は、ユークリッドの互除法というアルゴリズムが有名だ。ところでPythonにはfractionsモジュールの中に最大公約数を求める関数が組み込まれている。

from fractions import gcd
a, b = 12, 35
print gcd(a, b)
# => 1
# よって互いに素

さて、ABC-tripleを定義するとようやくABC予想が何であるか分かる。

予想: abc-triple (a, b, c) の全てが c < rad(abc)^2 を満たす。

http://ja.wikipedia.org/wiki/ABC%E4%BA%88%E6%83%B3

ただしABC予想の正確な表現は、

任意のε>0において、c < rad(abc)^(1+ε)を満たす

というものである。
ところで、ABC予想においてはquality of an ABC-triplesという指標があり、

q(a, b, c) = log(c) / log(rad(abc))

が 1 より大きいと、そのABC-tripleは小さな素数の積で成り立っていることを示す(このへんなるほどわからん状態)。

例題:ABC-tripleの探索

今回は、cが1000以下における、q(a, b, c)が 1 より大きいABC-tripleをすべて発見してみよう。
コードは以下のようになる。radやgcdは上で見てきた通りのものを使っているので省略。

from math import log
from operator import itemgetter

def q(a, b, c):
    return log(c) / log(rad(a * b * c))

def find_abc(limit=1000):
    abc_triples = []
    for a in range(1, limit - 2):
        for b in range(a + 1, limit - a):
            c = a + b
            if c < limit and gcd(a, b) == 1:
                qabc = q(a, b, c)
                if qabc > 1:
                    abc_triples.append((qabc, a, b, c))
    return abc_triples

if __name__ == "__main__":
    import pprint
    pp = pprint.PrettyPrinter()
    getfirst = itemgetter(0)

    limit = 1000
    abc_triples = find_abc(limit)
    abc_triples = sorted(abc_triples, key=getfirst, reverse=True)

    print "c < %d において" % limit
    print "q(abc) > 1 の ABC-tripleの数:", len(abc_triples)
    print "最もq(abc)の高いABC-triple:", abc_triples[0]
    print "すべてのABC-triple:"
    pp.pprint(abc_triples)

# 出力
# c < 1000 において
# q(abc) > 1 の ABC-tripleの数: 31
# 最もq(abc)の高いABC-triple: (1.4265653296335432, 3, 125, 128)
# すべてのABC-triple:
# [(1.4265653296335432, 3, 125, 128),
# (1.3175706029138485, 1, 512, 513),
# (1.3111012199262271, 1, 242, 243),
#          ...(省略)...

a, bの分布図は次のようになる(縦軸の範囲を横軸の半分に取っている)。

最もq(abc)の高いABC-tripleは、次のように素因数分解される。

q a b c
1.4266 3 5^3 2^7

このq(abc)が大きいと何が嬉しいのかはまだわかりません。