素数判定・列挙

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

  • $n$ が素数か判定する
    • 試し割り、ミラーラビン素数判定、etc.
  • 素数を列挙する
    • $n$ までの素数を列挙する
      • エラトステネスのふるい、etc.
    • 小さい方から無限に素数を列挙する
  • $n$ を素因数分解する
    • 試し割り、ポラード=ロー法、etc.

素数判定

$n$ が素数か判定する

試し割り

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

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

高速化その1 奇数に絞る

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

実装例

高速化その2 素数に絞る

後述のエラトステネスのふるい等で $\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$ でフェルマーテストを通過(素数じゃね?と判定)されてしまう数も存在する。

Miller-Rabin素数判定法

フェルマーテストを改良した確率的アルゴリズム。
後述するように、$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個でもよい

実装例

nまでの素数を列挙する

エラトステネスのふるい

はじめに $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}$ なので、オーダー記法では省略される。

実装例

省メモリ化の工夫

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

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

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

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

アトキンのふるい

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

上限 algorithm time[s]
$10^6$ 無限ジェネレータ 0.43433
エラトステネス 0.12945
エラトステネス(NumPy) 0.11291
$10^8$ 無限ジェネレータ 65.4283
エラトステネス 14.7402
エラトステネス(NumPy) 9.8643

素因数分解

試し割り

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

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

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

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

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

実装例

Pollard's rho algorithm

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

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

  • $n$ が同じ桁数の2種類の素数の積であるとき、$O(n^{1/4} polylog(n))$ ステップ実行することで発見確率が約50%となる
programming_algorithm/number_theory/prime_judge.txt · 最終更新: 2024/01/19 by ikatakos
CC Attribution 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0