素数判定・列挙

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

素数判定

nが素数か判定する

試し割り

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

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

高速化その1

はじめに2で試し割りする。

以降は$\sqrt{n}$までの3以上の奇数で割る。偶数で割り切れるんなら、その前に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

高速化その2

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

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

Miller-Rabin素数判定法

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

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

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

フェルマーの小定理より、任意の素数$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,19,23
  • $2^{64}$までなら2,3,5,7,11,13,17,19,23,29,31,37

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

# 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^48までの決定的アルゴリズムとして実装しているので、ランダム要素は無い
# ランダムとして用いる場合は、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)
    for i in check_list:
        if i >= n:
            break
        if not suspect(i, d, n):
            return False
    return True

nまでの素数を列挙する

エラトステネスのふるい

エラトステネスの篩

はじめに$n$までの全ての数を素数候補に入れておき、素数だとわかった数の倍数を候補から外していく。

  • メリット
    • nまでの素数を一気に調べられる
    • 1000万くらいまでなら十分実用的
  • デメリット
    • 1億とかになるとちょっと時間的にきつい
    • 要素数nの配列を作るので、nが大きくなるとメモリ的にもきつい

省メモリ化の工夫

「30」周期で考える。これは2,3,5の最大公約数。すると、この3数のいずれかで割りきれる数は素数でないため、情報を持っておく必要が無い。

3数のいずれでも割りきれない数は、30の倍数に$(1,7,11,13,17,19,23,29)$のいずれかを足した数しか無い。よって、30個→8個に情報を減らせる。

アトキンのふるい

Sieve of Atkin

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

未調査。

無限に素数を列挙する

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

streamとfilter

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

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

  • 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に登録されていれば、合成数。逆に無ければ素数。

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

調べている数調べる前のhの中身説明
3{}hは空(hのkeyに3は無い)なので3は素数
次の3の倍数である9を登録(奇数しか調べないので6は飛ばす。以降も奇倍数のみ登録)
5{9: 3}hのkeyに5は無いので5は素数。次の5の倍数である15を登録
7{9: 3, 15: 5}hのkeyに7は無いので7は素数。次の7の倍数である21を登録
9{9: 3, 15: 5, 21: 7}keyに9があるので9は合成数
h[9]=3の倍数であったことがわかる
次の3の倍数である15を登録しようとするが、登録済みなのでスキップ
次の3の倍数である21を登録しようとするが、登録済みなのでスキップ
次の3の倍数である27を登録
9をkeyから削除
11{27: 3, 15: 5, 21: 7}11は素数
13{27: 3, 15: 5, 21: 7, 33: 11}13は素数
15{27: 3, 15: 5, 21: 7, 33: 11, 39: 13}h[15]=5 があるので15は合成数
次の5の倍数である25を登録
15をkeyから削除
17{27: 3, 25: 5, 21: 7, 33: 11, 39: 13}17は素数
19{27: 3, 25: 5, 21: 7, 33: 11, 39: 13, 51: 17}19は素数

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

programming_algorithm/number_theory/prime_judge.txt · 最終更新: 2017/11/20 by ikatakos
CC Attribution 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0