素数判定
(math/is-prime.hpp)

View this file on GitHub
Last update: 2023-10-01 12:06:18+09:00
Include: #include "math/is-prime.hpp"

Miller–Rabin 素数判定法

$p$ を奇素数とする。

$1$ の平方根

$p$ を法とした剰余において、$1$ の平方根として $-1,1$ が存在する。逆に、$-1,1$ 以外に $1$ の平方根が存在しないこと背理法により示す。

$1$ の平方根を $x (\neq 1,-1) \pmod p$ とすると、

$x^2 \equiv 1 \pmod p$

$(x-1)(x+1) \equiv 0 \pmod p$

$p$ は素数なので、$x+1,x-1$ のどちらか一方は必ず $p$ で割り切れる。しかし、$x \neq 1,-1 \pmod p$ より $x+1,x-1$ は $p$ で割り切れない。よって矛盾が生じる。

以上より $\mathrm{mod}\space p$ において $1$ の平方根は $1, -1$ のみである。

強擬素数(strong pseudoprime)

$n$ を $3$ 以上の奇数とすると、奇数 $d$、正整数 $s$ を用いて $n-1 = d \cdot 2^s$ と表せる。

$n$ が素数ならば、全ての $a \in {0,1,2,…,n-1}$ について:

$a^d \equiv 1 \pmod n$
ある $0 \leq r \leq s-1$ に対して、 $a^{d \cdot 2^r} \equiv -1 \pmod n$

のいずれかが成立する。… (a)

証明

フェルマーの小定理より、

$a^{n-1} \equiv a^{d \cdot 2^{s}} \equiv 1 \pmod n$

よって、$a^{n-1} \pmod n$ の平方根は $-1$ または $1 \pmod n$ である。

平方根が $-1$ の時は 2番目の式が成立する。$1$ だった時はまた平方根を考える。このように平方根を取り続け、1度も $-1$ にならず $r=0$ となった時を考える。

この時 $a^d \equiv 1 \pmod n$ が成り立ち、1番目の式が成立する。

以上より示された。

(a) の対偶を取ると、「ある $a \in {0,1,2,…,n-1}$ において、$ a^d \not\equiv 1 \pmod n$ かつ、任意の $0 \leq r \leq s-1$ に対して $a^{d \cdot 2^r} \not\equiv -1$ が成り立つ」となる。

このような $a$ を見つけた場合、$n$ は合成数だと判定出来る。このような $a$ を $n$ の合成性の証拠という。証拠にならない $a$ を「強い嘘つき(strong liar)」といい、$n$ を底 $a$ についての「強擬素数(strong pseudoprime)」という。

正しく素数判定出来る確率

この方法では $n$ が素数の時は必ず素数であると判定出来るが、$n$ が合成数の時に誤って素数だと判定してしまうことがある。

任意の奇数の合成数 $n$ において、少なくとも $\frac{3}{4}$ が合成数の証拠になることが分かっている。つまり、各 $a$ について誤って素数だと判定してしまう確率は $\frac{1}{4}$ である。

よって、$a$ を独立にランダムに $k$ 回選んで強擬素数かどうかの判定を行った場合、誤って素数だと判定してしまう確率は最悪でも $4^{-k}$ である。

なお、$n < 2^{32}$ の時、$a=2,7,61$を調べれば必ず素数判定が出来ることが分かっている。

また、$n < 2^{64}$ の時、$a=2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37$ を調べれば十分であることが分かっている。

以上より、競プロの範囲ではほとんどの場合決定的なアルゴリズムなので安心して利用して良い。

Depends on

math/pow_mod.hpp

Verified with

Code

#pragma once

#include <array>

#include "pow_mod.hpp"

/*
ref: Fast Primality Testing for Integers That Fit into a Machine Word
Michal Foriˇsek and Jakub Janˇcina
*/

constexpr bool is_prime(int n) {
    if(n <= 1) return false;
    if(n == 2 or n == 7 or n == 61) return true;
    if((n & 1) == 0) return false;
    long long d = n - 1;
    while((d & 1) == 0) d >>= 1;
    constexpr std::array<int, 3> bases = {2, 7, 61};
    for(int a : bases) {
        long long t = d;
        long long y = pow_mod(a, t, n);
        while(t != n - 1 && y != 1 && y != n - 1) {
            (y *= y) %= n;
            t <<= 1;
        }
        if(y != n - 1 && (t & 1) == 0) return false;
    }
    return true;
}

constexpr bool is_prime(long long n) {
    if(n <= 1) return false;
    if(n == 2) return true;
    if((n & 1) == 0) return false;
    long long d = n - 1;
    while((d & 1) == 0) d >>= 1;
    constexpr std::array<long long, 7> bases = {
        2, 325, 9375, 28178, 450775, 9780504, 1795265022
    };
    for(long long a : bases) {
        a %= n;
        if(a == 0) continue;
        long long t = d;
        long long y = pow_mod(a, t, n);
        while(t != n - 1 && y != 1 && y != n - 1) {
            y = (__int128_t)y * y % n;
            t <<= 1;
        }
        if(y != n - 1 && (t & 1) == 0) return false;
    }
    return true;
}

#line 2 "math/is-prime.hpp"

#include <array>

#line 2 "math/pow_mod.hpp"

constexpr long long pow_mod(long long x, long long k, long long m) {
    long long res = 1;
    long long mul = (x >= 0 ? x % m : x % m + m);
    while(k) {
        if(k & 1) res = (__int128_t)res * mul % m;
        mul = (__int128_t)mul * mul % m;
        k >>= 1;
    }
    return res;
}
#line 6 "math/is-prime.hpp"

/*
ref: Fast Primality Testing for Integers That Fit into a Machine Word
Michal Foriˇsek and Jakub Janˇcina
*/

constexpr bool is_prime(int n) {
    if(n <= 1) return false;
    if(n == 2 or n == 7 or n == 61) return true;
    if((n & 1) == 0) return false;
    long long d = n - 1;
    while((d & 1) == 0) d >>= 1;
    constexpr std::array<int, 3> bases = {2, 7, 61};
    for(int a : bases) {
        long long t = d;
        long long y = pow_mod(a, t, n);
        while(t != n - 1 && y != 1 && y != n - 1) {
            (y *= y) %= n;
            t <<= 1;
        }
        if(y != n - 1 && (t & 1) == 0) return false;
    }
    return true;
}

constexpr bool is_prime(long long n) {
    if(n <= 1) return false;
    if(n == 2) return true;
    if((n & 1) == 0) return false;
    long long d = n - 1;
    while((d & 1) == 0) d >>= 1;
    constexpr std::array<long long, 7> bases = {
        2, 325, 9375, 28178, 450775, 9780504, 1795265022
    };
    for(long long a : bases) {
        a %= n;
        if(a == 0) continue;
        long long t = d;
        long long y = pow_mod(a, t, n);
        while(t != n - 1 && y != 1 && y != n - 1) {
            y = (__int128_t)y * y % n;
            t <<= 1;
        }
        if(y != n - 1 && (t & 1) == 0) return false;
    }
    return true;
}

素数判定 (math/is-prime.hpp)