容赦のない夏が来ました。毎日猛暑が続いていますが、みなさんいかがお過ごしでしょうか。

私はあまりに暑くて外に出たくないので、代わりに読書を捗らせることにしました。最近なかなか数学をやる時間がとれていなかったので、2年前にTeX界隈でいつもお世話になっている山本さん（@munepixyz）が組版を担当しているという理由で購入した整数論の本をぼちぼち読んでいくことにしました。

発見・予想を積み重ねる ―それが整数論

安福 悠

オーム社

Amazonで詳細を見る

私の専門は自然言語処理で、暗号学のように特に整数論が密接に関わる分野でもないですが（むしろ機械学習が幅を効かせている分野なので線形代数や解析学の方が大事かも）、道楽として数学をやるには整数論はなかなかロマンがあって良いなぁと昔から思っていました。あと、強いて言えば稀にCTFをやることがあるので（これも道楽です。筆者は全然強くありません）整数論の知識があるとチョット得することがあるかもなぁという気持ちも若干あります。

ただ読み進めるだけでは面白くないので、本書に登場する主張や証明を、気の向くままにコードを書くことで検証しながら読んでいくことにしました。言語は何でもよかったのですが、書捨てコードをサクッと書けてライブラリも充実していそうということで、Pythonを使うことにします。特にPythonに強いこだわりがあるわけではないので、適宜別の言語を使用する可能性もあります。

ナイーブな素数判定（試し割り）

私も伊達にOSS開発者をやっているわけではないので、ナイーブな素数判定ぐらい寝てても書けると思っていたのですが……意外と境界部分をきちんと処理する必要があって、最初に書いたコードをとりあえず動かしたところ“2 is NOT prime”と出力されるような始末でした。

具体的に「ナイーブな素数判定」というのは、ある整数$n$が素数であるかどうかを判定するのに

2以上$n$未満の整数について、順番に$n$を割り切れるかどうか確認する

という方法です（この愚直な方法は「試し割り」とも呼ばれるようです）。

「ナイーブな」と言っているので、上のアイデアをそのまま実装してもいいのですが、さすがにもう少し簡単にできる高速化の工夫をすることにします。そこで、上記の代わりに

3以上$\sqrt{n}$未満の奇数について、順番に$n$を割り切れるかどうか確認する

という方針を採用することにします。

まず、明らかに2以外の素数はすべて奇数です。そこで、2だけを例外扱いとして他は$n$が偶数である限りは「素数でない」判定を下すこととして、あとは奇数の約数を持たないかどうかだけを確認していきます（奇数は偶数の約数を持ちません）。これにより、ループの回数を半分に減らせます。

また整数$a$が$n$の約数のとき、$a\times b = n$なる$n$の約数$b$が存在します。このとき$a, b$のいずれかは$\sqrt{n}$以下なので¹、$\sqrt{n}$以下に約数がなければ、それ以降には約数は存在し得ないことがわかります。よって、$\sqrt{n}$未満の整数についてのみ調べれば十分です。

以上の方針をPythonコードに落としたものが次です：

from math import sqrt, floor

def is_prime(n):
    # prime is larger than 1
    if n < 2:
        return False

    # 2 is prime
    if n == 2:
        return True

    # can n be devided by 2?
    if n % 2 == 0:
        return False

    # can n be devided by an odd integer lower than sqrt(n)?
    for i in range(2, floor(sqrt(n)), 2):
        if n % (i + 1) == 0:
            return False

    # if not, n is a prime!
    return True

定数倍の高速化をしたところで時間計算量は$O(\sqrt{n})$なのですが、こんな実装でもインターネット上で利用可能な素数判定Webアプリケーションの中でも比較的大きな自然数を扱えるものが対応している16桁よりも大きな自然数（17桁）について、手元のMacBook Airで30秒ほどで結果を得ることができました。

>>> from prime import *
>>> is_prime(99999999999999997)
True

ところで、巷ではしばしば「Pythonのforループは遅い！」「代わりにリスト内包表記を使うと早い！」と言われていますが、これはいつでもどこでも真なわけではありません。

上記の素数判定コードは、forループの代わりにジェネレータ内包表記を用いて次のように書くこともできます²：

def is_prime2(n):
    # prime is larger than 1
    if n < 2:
        return False

    # 2 is prime
    if n == 2:
        return True

    # can n be devided by 2?
    if n % 2 == 0:
        return False

    # if n can be devided by an integer lower than sqrt(n), n is NOT prime;
    # otherwise n is a prime!
    return not 0 in (n % (i + 1) for i in range(2, floor(sqrt(n)), 2))

これでis_prime()とis_prime2()という2つの等価な関数が用意できたわけですが、それぞれについてIPythonを用いて実行時間（の平均）を計測してみると、次のようになりました。

%timeit is_prime(99999999999999997)
26 s ± 111 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit is_prime2(99999999999999997)
33.2 s ± 182 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

結果を見ると、内包表記の方が25%ほど遅いことがわかります。Pythonで「forループが遅くて内包表記だと早い」というのは、あくまで

ls = []
for i in range(n):
    ls.append(i)

のような場合の話なのであって（list.append()のオーバーヘッドが大きい）、Pythonではいつでもどこでも内包表記が正義というわけではないことがわかります。

エラトステネスの篩

素数判定のついでに、件の整数論本には載っていませんが、有名な「エラトステネスの篩」も実装してみました。先のis_prime()関数は自然数$n$を与えたときに$n$が素数かどうかを判定するものでしたが、エラトステネスの篩は「$n$以下の素数をすべて列挙する」効率的な手法です。

大層な名前が付いていますが、実際の手順としては

1. 2以上$n$以下の自然数を並べる
2. 印が付いておらず、かつまだ$p$として採用したことのない数のうち最小のものを$p$とおく 3.$p$以外の$p$の倍数に印をつける
3. 手順2.–3.を繰り返し、$p$が$\sqrt{n}$を超えた時点で操作を終了

という小学生の頃に誰もがやったことがあるだろう素数の発見方法です（小学生の頃は、終了条件はハッキリしなかったかもしれませんが）。

この操作を、なるべく素直にPythonで実装したのが以下です。マーク情報を保持するリストmarkedと別に整数列を保持するのは非効率なので、markedのインデックスがその役割を果たせるように調整しています。

from math import sqrt, floor


def sieve_of_eratosthenes(n):
    # initialize the mark list
    marked = [False] * (n + 1)

    # prime is larger than 1
    marked[0:2] = True, True

    # iterate from 2 to sqrt(n)
    for i in range(2, floor(sqrt(n)) + 1):
        # if i is marked, i is NOT prime: skip
        if marked[i]:
            continue

        # if i is unmarked, i is prime: mark the multiples
        for j in range(i * 2, n + 1, i):
            marked[j] = True

    # unmarked numbers are primes
    for i in range(n + 1):
        if not marked[i]:
            yield i

このコードの漸近計算量ですが、外側のループが$n$のオーダーなのは明らかです。内側のループはちょっと難しいですが、$\sqrt{n}$までの素数の逆数和のオーダーで、具体的には$\log\log n$になるようです³。よってコード全体の時間計算量は$O(n\log\log n)$、空間計算量はmarkedの大きさから明らかに$O(n)$となるでしょう。