素数判定アルゴリズムについて

この記事はCCS Advent Calendar 2020の11日目の記事として作成されました。

前の記事：
note.com

nullです。今回は素数判定や素因数分解を行うアルゴリズムについて色々書きます。

整数 $N$ の素数判定

最も初歩的なものとしては、 $2$ 以上 $N-1$ 以下の整数が $N$ を割り切るかどうか調べる、という方法が考えられます。
例えばC++での実装は以下のようになります。

#include <iostream>

int main()
{
  long long n;
  std::cin >> n;
  if (n <= 1)
  {
    std::cout << n << " is not prime.\n";
    return 0;
  }
  bool is_prime = true;
  for (long long i = 2; i <= n - 1; i++)
  {
    if (n % i == 0)
    {
      is_prime = false;
      break;
    }
  }
  if (is_prime)
    std::cout << n << " is prime.\n";
  else
    std::cout << n << " is not prime.\n";
}

このコードは最大で $N-2$ 個の候補について探索を行うため、（時間）計算量 $O(N)$ のコードといえます。
計算量とは、「あるアルゴリズムの効率がどれだけ優れているか」を表す指標です。この場合であれば探索の回数が大体 $N$ に比例することを意味しています。

C++であれば1秒に $10^9$ 回程度のループを回せるので、このコードでは $N \leq 10^9$ 程度であれば現実的な時間で素数判定が出来ます。

しかし $N$ がこれより大きい場合には途方もない時間がかかってしまいます。もっと効率のよい方法はないでしょうか。

$N \geq 2$ が合成数であることは、 $N = A \times B \ (A,B \geq 2)$ と2つの2以上の整数の積に分解できることと同値です。
ここで重要な考察として、 ${\rm min}(A,B) \leq \sqrt{N}$ が成り立ちます。
よって実は、 $2$ 以上 $\lfloor \sqrt{N} \rfloor$ 以下の整数が $N$ を割り切るかどうか調べるだけで十分であるとわかります。

#include <iostream>

int main()
{
  long long n;
  std::cin >> n;
  if (n <= 1)
  {
    std::cout << n << " is not prime.\n";
    return 0;
  }
  bool is_prime = true;
  for (long long i = 2; i * i <= n; i++)
  {
    if (n % i == 0)
    {
      is_prime = false;
      break;
    }
  }
  if (is_prime)
    std::cout << n << " is prime.\n";
  else
    std::cout << n << " is not prime.\n";
}

このコードの計算量は $O(\sqrt{N})$ であるため、 $N \leq 10^{18}$ 程度であれば判定できます。 $N$ の上限が先程の $10^9 = 1000000000$ 倍になりました！
このように計算量のよいアルゴリズムを考えることで、規模の大きい問題を現実的な時間で解くことが出来ます。

整数 $N_1, N_2, \dots, N_k$ の素数判定

今度は複数の整数の素数判定を行う場合を考えます。もちろん一つ一つの整数に対して先程のアルゴリズムを適用してもよいのですが、このような場合には、ある種の「前計算」を行うことで全体での判定効率がアップします。

${\rm max} N_i = M$ として、 $M$ 以下の素数を列挙したとしましょう。その場合、列挙した素数の中に $N_i$ が含まれるかどうかを調べるだけでよいので、各判定を $O(1)$ で行えます。魅力的ですね。

問題はどうやって素数を列挙するかですが、これを上手く実現するアルゴリズムとして、エラトステネスの篩が知られています。
$i = 2,3,\dots,M$ の順に、もし $i$ が素数（＝まだ消されていない）ならば $i$ を素数と確定した上で $2i,3i,\dots$ を消し、素数でない（＝既に消されている）ならば何もしない、という操作をすることで素数を列挙することが出来ます。

このアルゴリズムの計算量は $O(M \ {\rm log}{\rm log} M)$ であることが知られています。*1
このように、前計算を行うことで各判定を効率よく行うことが出来ます。

#include <iostream>

const int MAX_N = 1e7;
bool is_prime[MAX_N + 1];

void make_primes()
{
  std::fill(is_prime, is_prime + MAX_N + 1, true);
  for (int i = 2; i <= MAX_N; i++)
  {
    // 既に消されている場合はスルー
    if (!is_prime[i])
      continue;
    // i の倍数を消す
    for (int j = i + i; j <= MAX_N; j += i)
      is_prime[j] = false;
  }
}

int main()
{
  make_primes(); // 素数を前計算
  int n;
  while (std::cin >> n)
  {
    if (n >= 2 and is_prime[n])
      std::cout << n << " is prime.\n";
    else
      std::cout << n << " is not prime.\n";
  }
}

$N$ の素因数分解

一般に、巨大な数の素因数分解は非常に困難な問題とされています。
しかし、そこそこの大きさの数であれば、現実的な時間で素因数分解を行うことが可能です。

エラトステネスの篩を拡張してみましょう。
数を消す代わりに、その数が何で割り切れるかを記録することにします。そうすると前計算の計算量は変わらないまま、素因数分解を行うことができます。
再帰的に、 $N$ が素数でなければ $N$ を先程記録した数で割って同様の操作を繰り返すことで、素因数分解が行えます。

$N$ の素因数の個数は最大でも ${\rm log}_2 N$ 個程度であるため、計算量は $O({\rm log} N)$ です。
ただし実装の際は連想配列std::mapを用いているため、素因数分解の計算量は $O({\rm log^2} N)$ となっています。（std::unordered_mapを使うとならし $O({\rm log} N)$ になります）

#include <iostream>
#include <map>

const int MAX_N = 1e7;
int can_div[MAX_N + 1] = {}; // iを割り切る2以上i未満の素数（なければ0）

void make_primes()
{
  for (int i = 2; i <= MAX_N; i++)
  {
    // 既に消されている場合はスルー
    if (can_div[i] != 0)
      continue;
    // i の倍数を消す
    for (int j = i + i; j <= MAX_N; j += i)
      can_div[j] = i;
  }
}

// nの素因数分解
std::map<int, int> factorization(int n)
{
  if (can_div[n] == 0)
  {
    std::map<int, int> mp;
    mp[n] = 1;
    return mp;
  }
  std::map<int, int> mp = factorization(n / can_div[n]);
  mp[can_div[n]]++;
  return mp;
}

int main()
{
  make_primes(); // 素数を前計算
  int n;
  while (std::cin >> n)
  {
    std::cout << n << " = ";
    std::map<int, int> mp = factorization(n);
    for (auto itr = mp.begin(); itr != mp.end(); itr++)
    {
      auto [p, e] = *itr;
      std::cout << p << "^" << e;
      if (next(itr) != mp.end())
        std::cout << " * ";
      else
        std::cout << "\n";
    }
  }
}

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

世の中には、「100%ではないが、十分な確率で正しい解が求まる」アルゴリズムがあります。これを確率的アルゴリズムといいます。
素数判定を行う確率的アルゴリズムのうち代表的なものは2つあり、フェルマーテストとミラーラビン素数判定法です。

フェルマーテストは、フェルマーの小定理を利用した判定法です。フェルマーの小定理は、 $p$ が素数かつ $a$ と $p$ が互いに素のときに $a^{p-1} \equiv 1 \ ({\rm mod} \ p)$ が成り立つというものです。

$N$ をフェルマーテストで素数判定する流れは次のようになります。

まず、 $2$ 以上 $N-1$ 以下の整数 $a$ をランダムに選びます。
もし $a$ と $N$ が互いに素でないならば、 $a$ と $N$ は $2$ 以上 $N-1$ 以下の公約数を持つので、その時点で $N$ は素数ではありません。
互いに素な場合は、もし $N$ が素数ならフェルマーの小定理より $a^{N-1} \equiv 1 \ ({\rm mod} \ N)$ が成り立つはずです。よってこれが成り立たなければ素数ではありません。

$N$ が素数かもしれない場合は、再び整数 $a$ をランダムに選びます。これを十分な回数繰り返します。

このテストはあくまでも $N$ が素数であるための必要条件を調べているだけにすぎないので、正しい判定ができるとは限りません。しかし実際にはそれなりの確率で素数判定ができることが知られています。*2一方でカーマイケル数と呼ばれる合成数に対しては、どのような $a$ を選んでも素数と判定されてしまうという問題点があります。

ミラーラビン素数判定法は、フェルマーテストを改善した、より高い確率で正しい素数判定を行える判定法です。フェルマーテストと違ってカーマイケル数に対しても高確率で正しい判定を行うことができます。ここでは名前を紹介する程度にとどめておきます。（自分もよくわかってないので）

おわりに

素数に関するアルゴリズムの研究は古くからなされているだけあって、とても奥が深いです。現在も盛んに研究がなされています。興味を持った方は是非調べてみてください。
また、計算量のよいアルゴリズムを見つけるということの大切さが伝わればよいなと思っています。これはゲーム制作やアプリ開発など、どの分野でも有用だと思います。