Baillie-PSW素数判定法

概要

この記事ではbaillie2021strengtheningの記述を基にBaillie-PSW素数判定法について解説する．Baillie-PSW素数判定法とは，底$2$のミラー・ラビン素数判定法とリュカの強素数判定法を組み合わせた手法である．このアルゴリズムは高速であり，擬素数（素数と誤って判定される合成数）の割合が相当少ないことが予想されている．そのため，最近のソフトウェアにおける素数判定法の実装でよく用いられている．

この記事での前提知識は合同式のみである．また，最短で素数判定法にたどり着けるようにするため，証明や細かい話は全て省略している．気になる人は参考文献を参照のこと．

プログラムが載っていなかったり，脚注が増えていたりなど一部内容が異なるが，この記事の原稿となった PDF版 のノートも配布する．

準備

最大公約数

最大公約数はユークリッドの互除法を用いて，計算量$O(\log n)$で計算できる．Pythonでの実装例は次の通りである．

      def gcd(a, b):
    while b:
        a, b = b, a % b
    return a
    

平方剰余とルジャンドル記号，ヤコビ記号

たとえば，法$5$の下では
\begin{align*} 1^2\equiv1, \quad 2^2\equiv4, \quad 3^2\equiv4, \quad 4^2\equiv1, \pmod{5} \end{align*}
なので，$1$と$4$は平方剰余，$2$と$3$は平方非剰余である．

特に法が奇素数のとき，平方剰余は様々な性質を持つ．このことを表すために次の記法を導入する．

たとえば，$p=5$のとき，
\begin{align*} \left(\frac{1}{5}\right)=1, \quad\left(\frac{2}{5}\right)=-1, \quad\left(\frac{3}{5}\right)=-1, \quad\left(\frac{4}{5}\right)=1, \quad\left(\frac{5}{5}\right)=0, \end{align*}
となる．

ルジャンドル記号は法を奇素数に制限した．それを正の奇数に拡張したものがヤコビ記号である．

$n$が奇素数の場合はルジャンドル記号と一致する．ヤコビ記号は最速で計算量$O(M(n)\log n)$で計算できることが知られているbrent2010m．以下，$(a/n)$はヤコビ記号を表すものとする．

アルゴリズムは異なるが，Pythonを使った単純なアルゴリズムの実装例を紹介する．実際に動かしてみるとわかるが，大きい数でも十分高速に動作する．

Jacobi symbol - Rosetta Code

      def jacobi(a, n):
    if n <= 0:
        raise ValueError("'n' must be a positive integer.")
    if n % 2 == 0:
        raise ValueError("'n' must be odd.")
    a %= n
    result = 1
    while a != 0:
        while a % 2 == 0:
            a /= 2
            n_mod_8 = n % 8
            if n_mod_8 in (3, 5):
                result = -result
        a, n = n, a
        if a % 4 == 3 and n % 4 == 3:
            result = -result
        a %= n
    if n == 1:
        return result
    else:
        return 0
    

従来の素数判定法

フェルマーテスト

法が素数のとき，フェルマーの小定理（Fermat's little theorem）という有名な定理が成り立つ．

一般にフェルマーの小定理の逆は成立しない．実際，
\begin{equation*} 2^{561-1} \equiv 4^{561-1} \equiv \cdots \equiv 1 \pmod{561} \end{equation*}
であるが，$561=3 \times 11 \times 17$は合成数である．このように，$\gcd(a,n)=1$となる任意の整数$a$について，$a^{n-1}\equiv1\pmod{n}$を満たす合成数$n$をカーマイケル数（Carmichael's number）と呼び，無数に存在する．実際，$10^{21}$までのカーマイケル数は$20{,}138{,}200$個であるpinch2007carmichael．とはいえ，こういったフェルマーの小定理の逆を満たす合成数が稀であることを期待して，次のような確率的素数判定法が提案された．

上述の通り，フェルマーテストを潜り抜ける合成数$n$は無数に存在する．この合成数$n$を判定時の整数$a$を用いて，底$a$の擬素数（base-$\boldsymbol{a}$ pseudoprime）という．フェルマーテストは後述の確率的素数判定法と比べると精度が悪いので，実際に用いられることは滅多にない．

Pythonでの実装例は次の通りである．

      import math
import random

def is_prime(n, KMAX=100):
    if n == 1: return False
    if n == 2: return True
    for k in range(KMAX):
        a = random.randint(2, n - 1)
        if math.gcd(n, a) != 1:
            return False
        if pow(a, n - 1, n) != 1:
            return False
    return True
    

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

ミラー・ラビン素数判定法（Miller-Rabin primality test）は次の定理に基づいた確率的素数判定法である．

ミラー・ラビン素数判定法はこの定理の対偶を利用して，合成数を判定する．

フェルマーテストと違って「強（strong）」がついているのは，素数の確率がより高いことを意味する．同様に，素数と誤って判定されてしまう合成数$n$を，底$a$の強擬素数（base-$\boldsymbol{a}$ strong pseudoprime）という．精度と計算効率の良さから，実用上はミラー・ラビン素数判定法が採用されることが多い．ところが，最近ではこのミラー・ラビン素数判定法と後述のリュカの強素数判定法を組み合わせたBaillie-PSW素数判定法が使われることが増えている．

ミラー・ラビン素数判定法のPythonの実装例を紹介する．

素数判定いろいろ - フェルマーテスト・ミラーラビン素数判定法 #Python - Qiita

      import random

def is_prime(n):
    if n == 2: return True
    if n == 1 or n & 1 == 0: return False
    d = (n - 1) >> 1
    while d & 1 == 0:
        d >>= 1
    for k in range(100):
        a = random.randint(1, n - 1)
        t = d
        y = pow(a, t, n)
        while t != n - 1 and y != 1 and y != n - 1:
            y = (y * y) % n
            t <<= 1
        if y != n - 1 and t & 1 == 0:
            return False
    return True
    

Baillie-PSW素数判定法

リュカ数列とリュカの素数判定法

リュカ数列と素数に対して，次の定理が成り立つ．

リュカ数列のパラメータを$(D/p)=-1$となるように取ると，この定理は次のように書ける．
\begin{align} U_{p+1} &\equiv 0 \pmod{p}, \tag{3}\label{eq:Uprime}\\ V_{p+1} &\equiv 2Q \pmod{p}.\tag{4}\label{eq:Vprime} \end{align}
式(\ref{eq:Uprime})から，フェルマーテストと同様の論法でリュカの素数判定法を定義する．ちなみに，式(\ref{eq:dame1}), (\ref{eq:dame2})は合成数であっても成立するケースが多いので素数判定では使われない．

リュカの素数判定法で素数と誤って判定される合成数をリュカの擬素数（Lucas pseudoprime）という．また，$\{V_k\}$の式(\ref{eq:Vprime})についてのリュカの素数判定法で得られる確率的素数をリュカの$V$確率的素数（Lucas-V probable prime）と呼ぶ．これを使うと，計算量を悪化させることなくリュカの素数判定法の判定率を向上できる．しかし，現状はまだ普及していないので，この記事では紹介しない．

リュカの素数判定法におけるパラメータ$(P,Q)$の選び方の一例を紹介する．3つ目のステップが無くても，生成される擬素数は同じである．また，これは後述のリュカの強素数判定法においても同様である．

• $D$を$5, -7, 9, -11, 13, -15, \dots$という数列の中で$(D/n)=-1$を最初に満たす要素として選ぶ．
• $P=1$, $Q=(1-D)/4$とする．
• もし，$Q=-1$であれば，$P=Q=5$とする．

この方法を用いた場合，リュカの擬素数は
$$ 323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179, 10877,\dots $$

となる．また，$10^{15}$以下のリュカの擬素数は$2{,}402{,}549$個である．

フェルマーテストからミラー・ラビン素数判定法への強化と同じようにして，次の定理を用いてリュカの強確率的素数と強擬素数を定義する．

リュカの強素数判定法をすり抜ける合成数を，パラメータ$(P,Q)$に対するリュカの強擬素数（strong Lucas pseudoprime）と呼ぶ．リュカの素数判定法と同様の方法でパラメータを選択した場合，リュカの強擬素数は
$$ 5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, \dots $$
となる．また，$10^{15}$以下のリュカの強擬素数は$474{,}971$個であり，大幅に判定率が改善されていることがわかる．

Baillie-PSW素数判定法

Baillie-PSW素数判定法は，ミラー・ラビン素数判定法とリュカの強素数判定法を組み合わせた手法である．ミラー・ラビン素数判定法における強擬素数と，リュカの強擬素数は異なる数となる傾向がある．すなわち，両方の素数判定法をすり抜ける擬素数は極めて稀ということである．実際，$2^{64}$以下では擬素数が存在しないことが知られている．

実用上，ミラー・ラビン素数判定法の前にしきい値$k$以下の素数で試し割り法を行うことが多い．また，この素数判定法は非常に効率的なので，多くのソフトウェアではこの方法が採用されているalbrecht2018prime．

さきほど，Baillie-PSW素数判定法をすり抜ける合成数が極めて稀と述べたが，実際のところは未解決問題である．ヒューリスティックな解析によると，無数に擬素数は存在すると予想されているが不明であるpomerance1984there．もし，反例を見つけた場合は620ドルの賞金が手に入るようだ．ちなみに，文献によって書いている賞金額が違うので，実際にはもっと少ない可能性がある．baillie2021strengtheningでは620ドル，pomerance1984thereでは120ドルと書かれている．また wikipedia によると，620ドルは別の問題に対する賞金だという指摘もある．

Baillie-PSW素数判定法のPythonの実装例は以下のページを参照のこと．愚直に実装するとリュカ数列の計算量が大変なことになるが，ビット表示して効率的に計算できる．

Baillie-PSW/baillie_psw.py at main · armchaircaver/Baillie-PSW (github.com)