解説 | 素数判定 (64bit)

解説

まず、エラトステネスの篩で $\sqrt{2^{64}}=2^{32}=4294967296$ 以下の自然数に含まれる $203280221$ 個の素数を用意し、それを用いて試し割り法を用いる方法で間に合うかについて言及します。 1秒間に 10億回 = $10^9$ 回試し割りできるとして (1クエリあたり約0.2秒)、この問題では最大 $10000$ 個の整数が与えられるため、これでは実行時間が間に合いません。

手元の環境では、比較的高速な実装をしても $2^{32}$ 以下の素数を列挙するためのエラトステネスの篩に約0.7秒、高速な整除判定用の定数生成に約1.2秒、2億回の整除判定に 1ケースあたり約0.2～0.3秒程度の実行時間、高速な整除判定用の定数を保存するメモリ空間約3.2GiB 程度 ( $203280221 \times 16$ bytes ) + α を要しました。（実行時間・メモリ使用量ともに実行制限を超過）

そのため、より高速な素数判定法を用いる必要があります。

また、この問題で与えられる数は 64ビット符号なし整数型を正しく扱えているか、それが難しい言語でも工夫して実装できているか、をテストの一環に含めているため、問題のクエリごとの制約は $2^{64}$ 未満の非負整数としています。

この問題をクリアしたら、64ビット数の素因数分解にも挑戦してみてください。問題例: Library Checker: Factolize

想定解: Miller-Rabin素数判定法

$2^{64}$ 未満の範囲においては、 $\lbrace 2, 325, 9375, 28178, 450775, 9780504, 1795265022 \rbrace$ をそれぞれ底とする7回の Miller-Rabin素数テストを用いる事で、決定的に素数を判定できます。 $2^{32}$ 未満の範囲においては、 $\lbrace 2,7,61 \rbrace$ をそれぞれ底とする3回の Miller-Rabin素数テストを用いて素数判定できます。参考: Deterministic variants of the Miller-Rabin primality test, Created and maintained by Wojciech Izykowski.
また、与えられた数を適切なハッシュ関数で前処理する事で、 Miller-Rabin素数テストの回数を $2^{32}$ 未満の範囲では $1$ 回、 $2^{64}$ 未満の範囲では (1.6万～26万個の定数と引き換えに) $2$ ～ $3$ 回に抑える方法も知られています。参考: Forisek, Michal and Jakub Jancina. “Fast Primality Testing for Integers That Fit into a Machine Word.” Conference on Current Trends in Theory and Practice of Informatics (2015).
ミラー–ラビン素数判定法, ja.wikipedia.org

Writer解: C++ (GCC/Clang)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47#include <cstdio>
#include <cstdint>
// C++ : 累乗剰余の簡易実装、 GCC/Clang 128bit整数型を使用 (MSVCなどでは使用不可)
// (128bit整数除算は特に遅いので、可能ならモンゴメリ剰余乗算などで実装したほうが速い)
uint64_t modpow(uint64_t b, uint64_t p, const uint64_t n) {
    if(p == 0) { return (n > 1) ? 1 : 0; }
    while((p & 1) == 0) {
	    p >>= 1;
        b = (uint64_t)(((__uint128_t)b) * ((__uint128_t)b) % ((__uint128_t)n));
    }
    uint64_t r = b;
    for(p >>= 1; p != 0; p >>= 1) {
        b = (uint64_t)(((__uint128_t)b) * ((__uint128_t)b) % ((__uint128_t)n));
        if((p & 1) == 1) { r = (uint64_t)(((__uint128_t)r) * ((__uint128_t)b) % ((__uint128_t)n)); }
    }
    return r;
}
// C++ : nは非負整数とし、1以下の整数もしくは合成数ならFalse、素数ならTrueを返す
// ミラー・ラビン素数判定法の決定論的変種: https://miller-rabin.appspot.com/
bool miller_rabin(const uint64_t n) {
	if(n == 2) { return true; }
    if(n < 2 || (n & 1) == 0) { return false; }
    uint64_t n1 = n - 1, d = n - 1;
    int s = 0;
    while((d & 1) == 0) { d >>= 1; s += 1; }
    uint64_t bases[] = {2,325,9375,28178,450775,9780504,1795265022};
    for(uint64_t a : bases) {
        a %= n;
        if(a == 0) { continue; }
        uint64_t t = modpow(a, d, n);
        if(t == 1) { continue; }
        for(int j = 1; t != n1; ++j) {
            if(j >= s) { return false; }
            t = modpow(t, 2, n);
		}
    }
    return true;
}
int main(int argc, char *argv[]) {
    int t;
    scanf("%d", &t);
    for(int i = 0; i < t; ++i) {
        unsigned long long a;
        scanf("%llu", &a);
        printf("%llu %d\n", a, miller_rabin(a) ? 1 : 0);
    }
}

Writer解: Rust

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48// Rust : 累乗剰余の簡易実装
// (128bit整数除算は特に遅いので、可能ならモンゴメリ剰余乗算などで実装したほうが速い)
fn modpow(mut b: u64, mut p: u64, n: u64) -> u64 {
    if p == 0 { return (n > 1) as u64; }
    while (p & 1) == 0 {
        p >>= 1;
        b = ((b as u128) * (b as u128) % (n as u128)) as u64;
    }
    let mut r = b;
    p >>= 1;
    while p != 0 {
        b = ((b as u128) * (b as u128) % (n as u128)) as u64;
        if (p & 1) == 1 {
            r = ((r as u128) * (b as u128) % (n as u128)) as u64;
        }
        p >>= 1
    }
    r
}
// Rust : nは非負整数とし、1以下の整数もしくは合成数ならfalse、素数ならtrueを返す
// ミラー・ラビン素数判定法の決定論的変種: https://miller-rabin.appspot.com/
fn miller_rabin(n: u64) -> bool {
    if n < 3 || (n & 1) == 0 { return n == 2; }
    let (n1, d, s) = (n - 1, (n - 1) >> (n - 1).trailing_zeros(), (n - 1).trailing_zeros());
    'bases: for &a in &[2,325,9375,28178,450775,9780504,1795265022] {
        let a = a % n;
        if a == 0 { continue; }
        let mut t = modpow(a, d, n);
        if t == 1 || t == n1 { continue; }
        for _ in 1..s {
            t = modpow(t, 2, n);
            if t == n1 { continue 'bases; }
        }
      return false;
  }
  true
}
fn main() {
    use std::io::prelude::*;
    let stdin = std::io::stdin();
    let mut input = stdin.lock().lines();
    let t = input.next().unwrap().unwrap().parse::<usize>().unwrap();
    for _ in 0..t {
        let a = input.next().unwrap().unwrap();
        println!("{} {}", a, miller_rabin(a.parse().unwrap()) as u32);
    }
}

Writer解: Python

Pythonでは TLE を防ぐため、 input() ではなく sys.stdin などを用いて入力の最適化をした方が良いようです。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30import sys
# Python : nは非負整数とし、1以下の整数もしくは合成数ならFalse、素数ならTrueを返す
# ミラー・ラビン素数判定法の決定論的変種: https://miller-rabin.appspot.com/
def miller_rabin(n: int) -> bool:
    if n <= 2 or (n & 1) == 0:
        return (n == 2)
    n1 = n - 1
    s = ((n1 & (-n1)) - 1).bit_length() # s : (n - 1) を2進数表記した時に最後から続く0の数 (trailing_zeros_count)
    d = n1 >> s
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        if a % n == 0:
            continue
        t = pow(a, d, n) # = pow(a, d) % n
        if t == 1 or t == n1:
            continue
        for _ in range(s - 1):
            t = pow(t, 2, n) # = pow(t, 2) % n
            if t == n1:
                break # else節を実行せず次の for a in ... ループへ
        else: # breakでループを抜けず、普通に for _ in range() ... ループが終了した時
            # for文に付くelse節の説明:
            # https://docs.python.org/ja/3/reference/compound_stmts.html#the-for-statement
            return False
    return True

t = int(sys.stdin.readline().rstrip())

for _ in range(t):
    a = int(sys.stdin.readline().rstrip())
    print(a, int(miller_rabin(a)))

Writer解: Ruby

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29# Ruby : nは非負整数とし、1以下の整数もしくは合成数ならfalse、素数ならtrueを返す
# ミラー・ラビン素数判定法の決定論的変種: https://miller-rabin.appspot.com/
def miller_rabin(n)
    return n == 2 if n <= 2 || n & 1 == 0
    n1 = n - 1
    s = ((n1 & -n1) - 1).bit_length # s : (n - 1) を2進数表記した時に最後から続く0の数 (trailing_zeros_count)
    d = n1 >> s
    [2,325,9375,28178,450775,9780504,1795265022].each do |a|
        a %= n
        next if a == 0
        t, f = a.pow(d, n), true
        next if t == 1 || t == n1
        (s - 1).times do
            t = t.pow(2, n)
            if t == n1
                f = false
                break
            end
        end
        return false if f
    end
    return true
end

t = gets.to_i
t.times do
    a = gets.to_i
    printf("%d %d\n", a, miller_rabin(a) ? 1 : 0)
end

別解: Baillie-PSW素数判定法

底 $\lbrace 2 \rbrace$ ${2}$ の Miller-Rabin素数テストと、強いリュカ擬素数(strong Lucas pseudoprime)テストを組み合わせた、 Baillie-PSW素数判定法という方法もあります。これは、2つの素数テストのいずれでも probably prime (おそらく素数) と判定されてしまう合成数が $2^{64}$ $2^{64}$ 未満の範囲では存在しない事を利用しています。
- OEIS:A001262 Strong pseudoprimes to base 2 : 底 $\lbrace 2 \rbrace$ の Miller-Rabin素数テストで probably prime (おそらく素数) と判定される合成数 (擬素数)
- OEIS:A217255 Strong Lucas pseudoprimes : 強いリュカ擬素数テストで probably prime (おそらく素数) と判定される合成数 (擬素数)
Baillie–PSW primality test, en.wikipedia.org
Strengthening the Baillie-PSW primality test Baillie-PSW素数判定法の判定漏れの可能性をさらに減らす強化

高速化: モンゴメリ剰余乗算、 Barrett reduction

$a\times b\mod n$ のような剰余乗算の計算を、加減算・乗算・シフト演算を主に使って行うモンゴメリ剰余乗算を使うと、 C や Rust などの高速な言語においては、低速な剰余算の実行回数を低減させることにより、 Miller-Rabin素数判定法や Baillie-PSW素数判定法をより速く計算することができます。Python などの比較的低速な言語においては、低速な剰余算の実行回数を減らせるメリットよりも、代わりに必要となる加減算・乗算・シフト演算の演算回数の増加によるデメリットの方が大きくなるため、逆効果になる場合があります。
法が $2^{32}$ 未満の範囲では、 Barrett reduction と呼ばれる方法も使われています。例: C++向けの AtCoder Library における、 dynamic modint の剰余乗算の実装
モンゴメリ乗算, ja.wikipedia.org
Barrett reduction, en.wikipedia.org

難易度評価

外部サイトへ移動します

素数判定 (64bit)

解説