ミラー・ラビン素数判定法
(library/math/miller-rabin.hpp)

View this file on GitHub
Last update: 2025-03-03 00:33:48+09:00
Include: #include "library/math/miller-rabin.hpp"

素数判定の確率的アルゴリズムの 1 つ。

ふつうの決定性アルゴリズムは $O(\sqrt{n})$ かかる。 AKS 素数判定法は $O((\log n)^{6})$ かかるのでまあ遅い。

このアルゴリズムは (たぶん) $O(k \log n)$ で済む。ここで $k$ は判定に使う基底の個数。 64 bit 整数だと高々 12 あれば決定的に判定できるらしい。

理論的に

(参考: 離散アルゴリズム特論 #13)

キーアイデア

$p$ が素数で $a$ が $1 \le a < p$ を満たす整数ならば、 $a^{p-1} \bmod p = 1$ である。 (フェルマーの小定理)

証明はまあいろんなところにあるので省略。

$p$ を素数として、 $a$ を $1 \le a < p$ かつ $a^2 \bmod p = 1$ なる整数とする。このとき、 $a = 1$ または $a = p-1$ である。

証明:

$a^2 \bmod p = 1$ より $(a-1)(a+1) \bmod p = 0$
よって $(a-1) \bmod p = 0$ または $(a+1) \bmod p = 0$
つまり $a \bmod p = 1$ または $a \bmod p = -1$
$1 \le a < p$ より $a = 1$ または $a = p-1$

アルゴリズム

まず、 $n$ を奇数とする (偶数なら $n=2$ を判定すればいいので)。そしたら $n-1 = u \cdot 2^k$ と表現する (ただし $u$ は奇数)。

以下の数列 $b_0, …, b_k$ を考える:

\[b_i = \begin{cases} a^u \bmod n & i = 0 \\ b_{i-1}^2 \bmod n & i > 0 \end{cases}\]

このとき、以下の 2 つが成り立つ。

$b_j = 1$ なる $j$ が存在しなければ (特に $b_k = a^{n-1} \ne 1$ であるから)、 $n$ は合成数である。 (フェルマーの小定理から)
$b_j = 1$ なる 最小の $j$ を $j_{\text{min}}$ とする。 $j_{\text{min}} > 0$ かつ $b_{j_{\text{min}}-1} \ne n-1$ ならば、 $n$ は合成数である。 (キーアイデア 2 より)

この条件のどちらかを満たす $a$ を $n$ の A-witness という。 ($n$ が合成数であることを示すための証拠 (witness) の意味)

実は、 $n$ が合成数であるとき、 $\lbrace 2, 3, …, n-1 \rbrace$ のうち $n$ の A-witness は $\dfrac{3}{4}$ 以上存在するらしい。よって、一様ランダムに $a \in \lbrace 2, 3, …, n-1 \rbrace$ を取り出せば、 $n$ が合成数であることを確率的に判定できる。これを $k$ 回繰り返して、すべてで素数と判定されれば「素数」と出力すれば、合成数を素数と間違える確率は $\dfrac{1}{4^k}$ 以下である。なお、素数を合成数と間違える確率は $0$ である。

基底を選ぶ

拡張リーマン予想 (?) を使うと、 $2$ から $\lfloor 2 \ln^2 n \rfloor$ までの $a$ を全部試せばいいらしい。 $n$ が小さい場合 (RSA とかで使われるものより小さいの意、例えば $2^{64}$ とか) だと、試す A-witness は少なくて済む。 Wikipedia に載っているものを見ると、例えば $n < 2047$ なら $a = 2$ だけでいいらしい。

Required by

ポラード・ロー素因数分解法 (library/math/pollards-rho.hpp)

Verified with

Code

#pragma once

#include <cassert>
#include <concepts>
#include <functional>
using namespace std;

namespace asalib {
  namespace math {
    template<integral T>
    constexpr bool miller_rabin_test(T n) {
      if (n < 2) [[unlikely]]
        return false;
      if (n % 2 == 0) return n == 2;
      T d = n - 1, s = 0;
      while (d % 2 == 0) d >>= 1, ++s;
      assert(d % 2 == 1);

      static_assert(sizeof(T) == 4 || sizeof(T) == 8, "T must be 32-bit or 64-bit integer type");
      constexpr bool is_64bit = sizeof(T) == 8;

      function<T(T, T, T)> modpow = [&](T a, T n, T mod) {
        T res = 1;
        while (n) {
          if constexpr (is_64bit) {
            if (n & 1) res = (__uint128_t) res * a % mod;
            a = (__uint128_t) a * a % mod;
          }
          else {
            if (n & 1) res = (unsigned long long) res * a % mod;
            a = (unsigned long long) a * a % mod;
          }
          n >>= 1;
        }
        return res;
      };

      function<bool(T)> test = [&](T a) {
        T x = modpow(a, d, n);
        if (x == 0) [[unlikely]]
          return false;
        if (x == 1) [[unlikely]]
          return true;
        for (T _ = 0; _ < s; ++_) {
          T y;
          if constexpr (is_64bit)
            y = (__uint128_t) x * x % n;
          else
            y = (unsigned long long) x * x % n;
          if (y == 1 && x != 1 && x != n - 1) return false;
          x = y;
        }
        return x == 1;
      };

      function<bool(vector<T>)> tests = [&](vector<T> basis) {
        for (T a : basis)
          if (!test(a)) return false;
        return true;
      };

      if (n < 2'047) return test(2);
      if (n < 1'373'653) return tests({ 2, 3 });
      if (n < 9'080'191) return tests({ 31, 73 });
      if (n < 25'326'001) return tests({ 2, 3, 5 });
      if (n < 3'215'031'751) return tests({ 2, 3, 5, 7 });
      if (n < 4'759'123'141) return tests({ 2, 7, 61 });
      if (n < 1'122'004'669'633) return tests({ 2, 13, 23, 1662803 });
      if (n < 2'152'302'898'747) return tests({ 2, 3, 5, 7, 11 });
      if (n < 3'474'749'660'383) return tests({ 2, 3, 5, 7, 11, 13 });
      if (n < 341'550'071'728'321) return tests({ 2, 3, 5, 7, 11, 13, 17 });
      if (n < 3'825'123'056'546'413'051) return tests({ 2, 3, 5, 7, 11, 13, 17, 19, 23 });
      if (ull(n) <= 18'446'744'073'709'551'615ULL) return tests({ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 });
      // 18'446'744'073'709'551'616 は 2^64 であり、これ以上の数は 64 ビット整数型では表現できない
      // TODO: 多倍長整数を使う場合は if を追加する
      // TODO: C++23 以降なら __builtin_unreachable() -> std::unreachable() に変更可能
      __builtin_unreachable();
    }
  }  // namespace math
}  // namespace asalib

#line 2 "library/math/miller-rabin.hpp"

#include <cassert>
#include <concepts>
#include <functional>
using namespace std;

namespace asalib {
  namespace math {
    template<integral T>
    constexpr bool miller_rabin_test(T n) {
      if (n < 2) [[unlikely]]
        return false;
      if (n % 2 == 0) return n == 2;
      T d = n - 1, s = 0;
      while (d % 2 == 0) d >>= 1, ++s;
      assert(d % 2 == 1);

      static_assert(sizeof(T) == 4 || sizeof(T) == 8, "T must be 32-bit or 64-bit integer type");
      constexpr bool is_64bit = sizeof(T) == 8;

      function<T(T, T, T)> modpow = [&](T a, T n, T mod) {
        T res = 1;
        while (n) {
          if constexpr (is_64bit) {
            if (n & 1) res = (__uint128_t) res * a % mod;
            a = (__uint128_t) a * a % mod;
          }
          else {
            if (n & 1) res = (unsigned long long) res * a % mod;
            a = (unsigned long long) a * a % mod;
          }
          n >>= 1;
        }
        return res;
      };

      function<bool(T)> test = [&](T a) {
        T x = modpow(a, d, n);
        if (x == 0) [[unlikely]]
          return false;
        if (x == 1) [[unlikely]]
          return true;
        for (T _ = 0; _ < s; ++_) {
          T y;
          if constexpr (is_64bit)
            y = (__uint128_t) x * x % n;
          else
            y = (unsigned long long) x * x % n;
          if (y == 1 && x != 1 && x != n - 1) return false;
          x = y;
        }
        return x == 1;
      };

      function<bool(vector<T>)> tests = [&](vector<T> basis) {
        for (T a : basis)
          if (!test(a)) return false;
        return true;
      };

      if (n < 2'047) return test(2);
      if (n < 1'373'653) return tests({ 2, 3 });
      if (n < 9'080'191) return tests({ 31, 73 });
      if (n < 25'326'001) return tests({ 2, 3, 5 });
      if (n < 3'215'031'751) return tests({ 2, 3, 5, 7 });
      if (n < 4'759'123'141) return tests({ 2, 7, 61 });
      if (n < 1'122'004'669'633) return tests({ 2, 13, 23, 1662803 });
      if (n < 2'152'302'898'747) return tests({ 2, 3, 5, 7, 11 });
      if (n < 3'474'749'660'383) return tests({ 2, 3, 5, 7, 11, 13 });
      if (n < 341'550'071'728'321) return tests({ 2, 3, 5, 7, 11, 13, 17 });
      if (n < 3'825'123'056'546'413'051) return tests({ 2, 3, 5, 7, 11, 13, 17, 19, 23 });
      if (ull(n) <= 18'446'744'073'709'551'615ULL) return tests({ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 });
      // 18'446'744'073'709'551'616 は 2^64 であり、これ以上の数は 64 ビット整数型では表現できない
      // TODO: 多倍長整数を使う場合は if を追加する
      // TODO: C++23 以降なら __builtin_unreachable() -> std::unreachable() に変更可能
      __builtin_unreachable();
    }
  }  // namespace math
}  // namespace asalib

ミラー・ラビン素数判定法 (library/math/miller-rabin.hpp)