素数判定・列挙

以下の内容についてまとめている。

• $n$ が素数か判定する
  • 試し割り、ミラーラビン素数判定、etc.
• 素数を列挙する
  • $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$で試したときに割り切れてるはず。無駄な試行を無くせる。

• フェルマーテスト - Wikipedia

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

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

フェルマーの小定理より、任意の素数 $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$ でフェルマーテストを通過（素数じゃね？と判定）されてしまう数も存在する。

• カーマイケル数

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

フェルマーテストを改良した確率的アルゴリズム。
後述するように、$n$ の上限を決めることで決定的アルゴリズムとすることもできる。

$n$ が素数かどうかを調べるにあたり、$n$ が偶数なら簡単なので、以下、奇数とする。

$p$ が素数の場合、$a^2 \equiv 1 \pmod{p}$ なる $a$ に対して、$a \equiv \pm 1$ となる。

• $(a-1)(a+1) \equiv 0$ と因数分解でき、$a-1,a+1$ のどちらかが0になるしかない

ここで、もし $p$ が素数で無く $p=qr$ とできる場合、 $(a-1)$ が $q$、$(a+1)$ が $r$ の倍数だったりしても等式は成立するので、必ずしも $a \equiv \pm 1$ とは限らなくなる。
これを利用する。

ミラーラビン素数判定法は以下の通りとなる。

• $n$ が2なら素数。それ以外の偶数なら合成数
• ある数 $a$ をいくつか決めて、以下を試す
  • $a^{n-1} \not\equiv 1$ の場合、$n$ は合成数
  • $d=(n-1)/2$ で初期化して、$d$ が偶数、かつ $a^d \equiv 1$ である間、$d←d/2$ とし続ける
    • $d$ が奇数になっても $a^d \equiv 1$ だった場合、$a$ を使ってはどちらとも断定できない
    • ある $d$ ではじめて1以外の数になった場合、
      • $a^d \equiv -1$ になった場合、$a$ を使ってはどちらとも断定できない
      • $a^d \not\equiv \pm 1$ になった場合、$n$ は合成数
• $a$ のうち1つでも合成数と判定されたら、合成数

以上で、フェルマーテストよりも、1つの $a$ からより多くの情報を得ることができる。

実装上は、$d$ を半分にしていく方法では累乗を求める部分で毎回 $O(\log{d})$ の計算が発生するが、 逆から調べる、つまり $d$ を2で割れるとこまで割ったところから始めて $d←2d$ としていくと、 前回の結果を2乗していけばいいので、若干計算量が減らせる。

$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までの素数）
  • または、$2, 325, 9375, 28178, 450775, 9780504, 1795265022$ の7個でもよい

• 参考:
  • Deterministic variants of the Miller-Rabin primality test. Miller-Rabin SPRP bases records
  • 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, 325, 9375, 28178, 450775, 9780504, 1795265022)
    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

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秒程度。

やっぱり基本となるのはこれ。

エラトステネスのふるいなどで、予め $\sqrt{n}$ までの素数を列挙しておく。

各素数 $p$ に対し、$n$ が割り切れる限り、素因数に $p$ を追加し、$n←n/p$ とする。

最後までやって、$n \neq 1$ なら、$n$ も素因数の1つとなる。

素数の列挙で $O(\sqrt{n} \log{\log{n}})$ かかるので、$n \ge 10^{15}$ くらいから数秒ではきつくなっていく。

実装例

def prime_factorize(x, primes=None):
    """ xを素因数分解、辞書型で返す """
    # primes: sqrt(x)までの素数リスト（昇順）
    #   省略時は奇数を試し割り（2だけ別処理）
    #   単発なら不要だが、複数のxで行う場合は事前計算した方がよい

    factors = {}
    if primes is None:
        if x & 1 == 0:
            b = (x & -x).bit_length() - 1
            if b > 0:
                factors[2] = b
                x >>= b
        primes = range(3, int(x ** 0.5) + 1, 2)
    for p in primes:
        if p * p > x:
            break
        if x % p == 0:
            factors[p] = 1
            x //= p
            while x % p == 0:
                factors[p] += 1
                x //= p
    if x != 1:
        factors[x] = 1
    return factors

• ポラード・ロー素因数分解法
• ポラード・ロー素因数分解法について - Qiita

1. ミラーラビン素数判定法 によって素数判定する
2. 合成数の場合、フロイドの循環検出法、または改良版のブレントの循環検出法により2つの約数に分解する
3. どの約数も素数判定されるまで、再帰的に繰り返す

ミラーラビンは確率的アルゴリズムだが、上限を決めると決定的にできる。
フロイドまたはブレントの循環検出法は、確率的アルゴリズムであり、以下のことが経験的に知られているが、正確な解析は未解決らしい。

• $n$ が同じ桁数の2種類の素数の積であるとき、$O(n^{1/4} polylog(n))$ ステップ実行することで発見確率が約50%となる