素数判定・列挙

• $n$ が素数か判定する
  • 試し割り、ミラーラビン素数判定、etc.
• $n$ までの素数を列挙する
  • エラトステネスのふるい、etc.
• 無限に素数を列挙する

素数判定

$2～\sqrt{n}$ までの整数でとにかく割ってみて、どれでも割り切れなければ素数。

愚直な方法だが、よほど $n$ が大きくない限り十分実用的。1000万とか億とかになってくると、ちょっと辛いので、高速な判定法を使う。

はじめに2で試し割りし、以降は $\sqrt{n}$ までの3以上の奇数で割る。
偶数で割り切れるのなら、その前に2で割り切れてるはず。単純に2倍高速。

# python
# 高速化その1で実装
def trial_division(n):
    if n == 2:
        return True
    if n < 2 or n % 2 == 0:
        return False
    for i in range(3, int(ceil(sqrt(n))) + 1, 2):
        if n % i == 0:
            return False
    return True

後述のエラトステネスのふるい等で $\sqrt{n}$ までの素数リストを作って、それで割ってみる。

その1と同様に、$p \times m$で割り切れるんなら、$p$で試したときに割り切れてるはず。無駄な試行を無くせる。

ミラー–ラビン素数判定法

確率的アルゴリズム。「素数じゃない」と判定された数は絶対素数じゃ無いが、「素数かも」と判定された数はごく稀に素数じゃない可能性がある。

そのようなデメリットがあれど、試し割りなどと比較してかなり高速に判定できる。

フェルマーの小定理より、任意の素数 $p$ について、$a$ と $p$ が互いに素であれば $a^{p-1} \equiv 1 \pmod{p}$ が成立する。
その対偶より、$a^{n-1} \not\equiv 1 \pmod{n}$ なら、nは素数じゃない、という発想が基本の部分。

だからといって「$a^{n-1} \equiv 1 \pmod{n}$ なら $n$ は素数」とは言えないが、 $a$ を変えて複数回試すことで「ここまでいろんな $a$ で成立するなら素数じゃね？」ということにする。 $a$ を $k$ 種類試した時の誤判定の確率は $4^{-k}$ とのことなので、50回も試せばまず間違い無いと言える。

$a$ の試す値を「ある範囲までなら正しい結果が得られることが確認されている要素」に固定することで、範囲を限定した決定的アルゴリズムとして用いることもできる。

• $2^{32}$ までなら2,7,61
• $2^{48}$ までなら2,3,5,7,11,13,17（17までの素数）
• $2^{64}$ までなら2,3,5,7,11,13,17,19,23,29,31,37（37までの素数）
• $4.75 \times 10^{9}$ までなら2,7,61
• $3.18 \times 10^{23}$ までなら37までの素数

目安として、上記を調べれば、正しいことが確認されているらしい。

• A014233 - OEIS
  • $n$ 番目の素数までで判定したときに、最初に誤判定される偽素数
• Miller–Rabin_primality_test

# python
# 補助関数
def suspect(a, t, n):
    x = pow(a, t, n)
    n1 = n - 1
    while t != n1 and x != 1 and x != n1:
        x = pow(x, 2, n)
        t <<= 1
    return t & 1 or x == n1

# メイン
# 2^64までの決定的アルゴリズムとして実装しているので、ランダム要素は無い
# ランダムとして用いる場合は、check_listにランダム抽出された数を採用し、20～50個程度試す
def miller_rabin(n):
    if n == 2:
        return True
    if n < 2 or n % 2 == 0:
        return False
    d = (n - 1) >> 1
    while d & 1 == 0:
        d >>= 1
    check_list = (2, 7, 61) if n < 2 ** 32 else (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)
    for i in check_list:
        if i >= n:
            break
        if not suspect(i, d, n):
            return False
    return True

• エラトステネスの篩

はじめに $n$ までの全ての数を素数候補に入れておき、素数だとわかった数の倍数を候補から外していく。
具体的な手順は、上記Wikipediaが詳しい。

• メリット
  • $n$ までの素数を一気に調べられる
  • 1000万くらいまでなら一瞬
• デメリット
  • それ以上の桁数となると、時間・メモリともに $n$ に比例してきつくなっていく

偶数は2以外素数じゃないの。当たり前だけどね。

だから、単純な枝刈りとして、2の倍数を除いた後は奇数のみチェックすることにすると参照すべき範囲を約半分に減らせる。

また、素数かどうか調べるのは、$\sqrt{n}$ までの範囲でよい。
もし $\sqrt{n}$ より大きい数が素数じゃないなら、既にそれ以下の数で割り切れ、候補から除かれている。
$\sqrt{n}$ まで調べた後に残っているものは、素数確定。

また、$p$ が素数とわかってその倍数を候補から外す際も、調べるのは $2p,3p,4p,...$ ではなく $p^2, p^2+2p, p^2+4p,...$ でよい。
$p^2$ 未満の $p$ の倍数は、既に $p$ 以下の素数をチェックした際に除かれている。
たとえば $7$ が素数とわかった場合、14,21,35などは既に2,3,5などをチェックした際に除かれているので、まだ除かれていない可能性があるのは $49$ からとなる。

計算量は、$O(n \log\log{n})$ となる。

• エラトステネスのふるいとその計算量 | 高校数学の美しい物語

見つかった素数の倍数を候補から除いていくのに全体で $\dfrac{n}{2}+\dfrac{n}{3}+\dfrac{n}{5}+...$ と、 $\sqrt{n}$ までの素数の逆数の和 $\times n$ の計算が必要になるが、これは上から $n\log\log{n}$ で押さえられるらしい。

実際には $\sqrt{n}$ までの要素が素数かチェックするのにも別途 $O(\sqrt{n})$ かかるが、$\sqrt{n} \lt n\log\log{n}$ なので、オーダー記法では省略される。

def sieve_eratosthenes(n):
    primes = [0, 1] * (n // 2 + 1)
    if n % 2 == 0:
        primes.pop()
    primes[1] = 0
    primes[2] = 1
    for p in range(3, int(n ** 0.5) + 1, 2):
        if primes[p]:
            for q in range(p * p, n + 1, 2 * p):
                primes[q] = 0
    return primes

NumPyを使うと以下の感じ（numpyモジュールの読み込みを除くと、10%ほど速い）

import numpy as np

def sieve_eratosthenes(n):
    primes = np.zeros(n + 1, dtype=np.bool)
    primes[2] = 1
    primes[3::2] = 1
    for p in range(3, int(n ** 0.5) + 1, 2):
        if primes[p]:
            primes[p * p::2 * p] = 0
    return primes

先ほどは2の倍数をチェックしないことで高速化を図ったが、そもそも primes 配列で情報を持たないようにすると、メモリも約半分になる。

また、2に加え3の倍数とわかっている数も飛ばすようにすると、6で割ってあまりが1か5となる数さえ調べればよく、チェックすべき数が元の1/3になる。

さらに5の倍数を除くと、30の倍数に $(1,7,11,13,17,19,23,29)$ のいずれかを足した数だけ調べればよい。約1/4に減らせる。

ただし、配列の何番目の項が実際は何を指すかの計算が生じるため、計算時間は単純にはこの通りには減らない。

Sieve of Atkin

かなーり複雑になったが、高速になったエラトステネスの改良版。

未調査。

呼ばれるたびに、次の素数を無限に返せるジェネレータ。

ただ、なんだかんだ仮にでも上限を決めて、エラトステネスの篩を適用した方が速いと思う。

実用性という意味では薄いが、まぁあまり大きな数が来ないとわかっているときに……（あるのか？）

pythonの、シーケンスとフィルターを使用する方法。

この辺で語られているのを発見。

• Python3で素数の無限ジェネレータ - Qiita

• itertools.countで無限に奇数を生成するジェネレータを作成
• 素数$p$が見つかると、以降「$p$で割り切れない」というフィルターを順次ジェネレータに重ねがける
  • よって、ジェネレータから取り出される値は必ず素数
• 記述こそ短いものの、あまり速くない
  • フィルターが何を行っているかというと、ジェネレータの各要素に対し、その数以下の全素数で試し割りしてる状態
  • 実際はsqrt(n)まででいいので、無駄が多い

def eratosthenes_stream():
    import itertools
    yield 2
    stream = itertools.count(3, 2)

    while True:
        prime = next(stream)
        stream = filter(lambda x, y=prime: x % y, stream)
        yield prime

• 汎用的で省メモリかつかなり速い素数無限列挙 - 名古屋で数学するプログラマ(仮)

こちらの方が圧倒的に高速。（上記サイトのrubyのソースをpythonにしただけ）

def eratosthenes_generator():
    yield 2
    n = 3
    h = {}
    while True:
        m = n
        if n in h:
            b = h[n]
            m += 2 * b
            while m in h:
                m += 2 * b
            h[m] = b
            del h[n]
        else:
            m += 2 * n
            while m in h:
                m += 2 * n
            h[m] = n
            yield n
        # print(n, h)
        n += 2

はじめに2を返した後は、奇数を順番に調べる。

特徴的なのは、dictを使用して情報管理している点。発見済みの各素数に対し「次に見つかる自分の倍数」をキー、自分を値として保持・更新していく。

よって、現在調べている数が既にdictのkeyに登録されていれば、合成数。逆に無ければ素数。

処理の順番は、以下のようになる。エラトステネスのふるいでの「次に見つかる倍数」のみを上手く管理していることがわかる。

1000万個の素数の列挙で、100秒程度。