Solovay-Strassenの素数判定法

与えられた数 $n$ が素数かどうかを誤り率 $2^{-k}$ 以下、時間 $O(k\log(n))$ で求める Solovay–Strassen 素数判定法を紹介します。ミラー・ラビン素数判定法と比べてやや単純で面白いと思います。平方剰余の相互法則がアルゴリズムに生かせるというのが熱い。

1 平方剰余/非剰余と $x^{(p-1)/2}$ の関係
2 ルジャンドル記号の二つの計算
3 ヤコビ記号は速い
4 アルゴリズム
5 正当性の証明
6 実装
7 参考文献

アルゴリズムのアイディアを理解するには、平方剰余/非剰余の性質を少し知っている必要があります。

平方剰余/非剰余と $x^{(p-1)/2}$ の関係

奇素数 $p$ に対して、$x^2=a \pmod p$ が解をもつとき $a$ を平方剰余、解を持たないとき $a$ を非平方剰余といいます。$(\mathbb{Z}/p\mathbb{Z})^\times$ は巡回群であり、任意の $a \in (\mathbb{Z}/p\mathbb{Z})^\times$ に対して $a=g^i$ なる $i$ が存在する $g$ （原始根）が存在します。つまり、$(\mathbb{Z}/p\mathbb{Z})^\times=\{g^0,g^1,\ldots,g^{p-2}\}$ です。原始根を用いて平方剰余は $g^{2i}$、非平方剰余は $g^{2i+1}$ と書け、$p$ が奇数より両者は同数存在します。

$x^2 = 1 \Leftrightarrow x=\pm1$ と $x^{p-1}=1$ より、$x^{(p-1)/2}=\pm1$ です。平方剰余は $g^{2i}$ の形で表せたことから、 $x^{(p-1)/2} = 1$ の解です。高々 $(p-1)/2$ 個しか解が存在しないことから、平方剰余で解は尽くされており、残りの平方非剰余は $x^{(p-1)/2} = -1$ の解です。

ルジャンドル記号の二つの計算

$a^{(p-1)/2}=\pm1$ と平方剰余/非剰余の関係を念頭に置いて、奇素数 $p$ と整数 $a$ に対して、ルジャンドル記号を$$\left(\frac{a}{p}\right) = \begin{cases} 1 &( a \equiv \eta^2 \pmod p \text{ なる } \eta \neq 0 \pmod p\text{ が存在})\\
-1 &( a \equiv \eta^2 \pmod p \text{ なる } \eta \text{ が存在しない})\\
0 & (a=0 \pmod p)
\end{cases}$$で定義します。平方剰余/非剰余が $g$ の冪の偶奇で決まったことを思い出すと、$\left(\frac{ab}{p}\right) = \left(\frac{a}{p}\right) \left(\frac{b}{p}\right)$ が成り立ちます。異なる奇素数 $p, q$ 同士のルジャンドル記号には平方剰余の相互法則$$\left(\frac{q}{p}\right) \left(\frac{p}{q}\right)= (-1)^{\frac{p-1}{2}\cdot\frac{q-1}{2}}$$が成り立つことが知られています。また $\left(\frac{2}{p}\right) = (-1)^{(p^2-1)/8}$ が成り立ちます。以上から、ルジャンドル記号の計算は $\left(\frac{a}{p}\right)=a^{(p-1)/2}$ を用いることのほかに、次のような方法があります。

$a=2b$ ならば、$\left(\frac{2b}{p}\right) = (-1)^{(p^2-1)/8}\left(\frac{b}{p}\right)$ とします。
$a$ が奇数ならば、平方剰余の相互法則 $\left(\frac{a}{p}\right) = (-1)^{\frac{p-1}{2}\cdot\frac{a-1}{2}}\left(\frac{p}{a}\right)$ を用いて、法を交換します。$\left(\frac{p}{a}\right) = \left(\frac{p \bmod a}{a}\right) $ より、法を交換すると底が必ず減らせます。法を交換したルジャンドル記号を再帰的に計算します。
$a$ が合成数ならば、素因数分解 $a=\prod p_i^{e_i}$ に対して、$\left(\frac{a}{p}\right) = \prod \left(\frac{p_i}{p}\right)^{e_i}$ です。各素因数のルジャンドル記号を再帰的に計算します。

平方剰余の相互法則により、$\left(\frac{a}{p}\right) = a^{(p-1)/2}$ とは全く異なった計算でルジャンドル記号を求められました。単にルジャンドル記号を求めるだけなら、平方剰余の相互法則を用いることは何の長所もありません（むしろ、因数分解で計算量が悪化）。しかし、$p$ の素数判定となれば、話は別です。「$ \left(\frac{a}{p}\right) $ を $ a^{(p-1) / 2} $ と平方剰余の相互法則の二通りで求めて、一致するか確かめる」というのが Solovay–Strassen の素数判定のアイディアです。

ヤコビ記号は速い

ルジャンドル記号を平方剰余の相互法則で求めるとき、素因数分解が計算量のボトルネックになりました。しかし、素因数分解は法を合成数に拡張したヤコビ記号で回避できます。ヤコビ記号は、任意の奇数 $n=\prod p_i^{e_i}$ と任意の整数 $a$ に対して、ルジャンドル記号の積 $\left(\frac{a}{n}\right) = \prod \left(\frac{a}{p_i}\right)^{e_i}$ で定義されます。ルジャンドル記号では、法が素数でなくとも、互いに素ならば平方剰余の相互法則が成り立ち、$$\left(\frac{m}{n}\right) \left(\frac{n}{m}\right)= (-1)^{\frac{n-1}{2}\cdot\frac{m-1}{2}}, \ \ \ \ (\gcd(n,m)=1)$$です。$n$ が合成数でも $\left(\frac{2}{n}\right)=(-1)^{(n^2-1)/8}$ は成り立ちます。$\gcd(n,m) \neq 1$ ならば $\left(\frac{n}{m}\right)=0$ であることを併せると、平方剰余の相互法則を繰り返し適用して法を交換することで、ヤコビ記号が $O(\log(n))$ で求まります！　ルジャンドル記号だけでは素因数分解が必要でしたが、ヤコビ記号では不要になり、計算量が改善しました。

アルゴリズム

アルゴリズムはシンプルです。誤り率 $2^{-k}$ 以下で $n$ が素数か判定したいとします。

$n$ が偶数ならば、$n=2$ のとき素数、さもなくば合成数です。以下、$n$ が奇数とします。
(i), (ii), (iii) を $k$ 回繰り返します。

ランダムな整数 $a \neq 0 \pmod n$ を取ってきます。
$\gcd(a,n) > 1$ ならば $n$ は合成数で確定です。
$\left(\frac{a}{n}\right)$ を平方剰余の相互法則、$a^{(n-1)/2}$ を繰り返し二乗法で求めます。$\left(\frac{a}{n}\right) \neq a^{(n-1)/2}$ ならば $n$ は合成数で確定です。

$n$ が合成数と確定していなければ、素数と判定します。

なかなか、シンプルです。合成数 $n$ を誤って素数と判定する確率が $2^{-k}$ 以下であることは次の節で証明します。

正当性の証明

$\left(\frac{a}{n}\right) \neq a^{(n-1)/2}$ ならば $n$ は合成数です。そこで、$a$ を合成数 $n$ に対するオイラー証人（Euler Witness）と呼びます。また、合成数 $n$ に対して $\left(\frac{a}{n}\right) = a^{(n-1)/2}$ となる $a$ をオイラー嘘つき（Euler Liar）と呼びます。Solovay–Strassen 素数判定法の誤り率が $2^{-k}$ 以下であることをこの節で示します。つまり、$\{1,2,\ldots,n-1\}$ のうち、オイラー証人が半分以上を占めることを示します。

定理：任意の合成数 $n$ に対して、オイラー証人は１つ以上存在します。

証明：
(i) $n$ が平方因子 $p^2$ ($p$ : 素数) を持つ場合：
$n$ と互いに素な数 $a$ に対して、定義より $\left(\frac{a}{n}\right)=\pm1$ です。従って、オイラー証人が存在しなければ任意の $a \in (\mathbb{Z}/n\mathbb{Z})^\times$ に対して $\left(\frac{a}{n}\right)^2=a^{n-1}=1$ です。$(\mathbb{Z}/p^2\mathbb{Z})^\times$ は要素数 $p(p-1)$ の巡回群なので、$n-1$ は $p(p-1)$ の倍数である必要があります。しかし、$n-1 \neq 0 \pmod p$ より矛盾しているので、オイラー証人が存在します。

(ii) ２つ以上の相異なる素数の積 $n=\prod p_i$ の場合：
中国剰余定理より、$i=1$ のとき $\left(\frac{a}{p_i}\right)=-1$、$i\neq 1$ のとき $\left(\frac{a}{p_i}\right)=1$ となる $a$ を取れます。このとき、$\left(\frac{a}{n}\right)=\prod \left(\frac{a}{p_i}\right)=-1$ です。同様に全ての $i$ で $\left(\frac{b}{p_i}\right)=1$ かつ $i \neq 1$ で $a = b \pmod {p_i}$ となる $b$ を取れます。このとき、$\left(\frac{b}{n}\right)=\prod \left(\frac{b}{p_i}\right)=1$ です。オイラー証人が存在しないと仮定します。$a^{(n-1)/2}=-1,b^{(n-1)/2}=1 \pmod n$ ですが、$a,b$ は $i \neq 1$ のとき $a^{(n-1)/2} = b^{(n-1)/2} \pmod {p_i}$ です。矛盾しているのでオイラー証人が存在します。

先の定理より、オイラー証人が一つは存在するので、次の定理より、オイラー証人が半分以上を占めます。

定理：任意の合成数 $n$ に対して、オイラー証人が一つでもあれば、$\{1,2,\ldots,n-1\}$ の半分以上がオイラー証人です。

証明：$x \neq 0 \pmod n$ とします。 $\gcd(x,n) > 1$ ならば、$n$ は合成数と確定するので、$x \not \in (\mathbb{Z}/n\mathbb{Z})^\times$ はオイラー証人です。オイラー嘘つき全体は $(\mathbb{Z}/n\mathbb{Z})^\times$ の部分群をなします。よって、ラグランジュの定理よりオイラー証人が存在すれば、オイラー証人は $(\mathbb{Z}/n\mathbb{Z})^\times$ の半分以上を占めます。

実装

短く書けてイイ感じ。


#include <random>
using namespace std;

long long pow(long long a, long long n, long long p) {
  return n != 0 ? pow(a*a%p, n/2, p) * (n%2 == 1 ? a : 1) % p : 1;
}

long long gcd(long long a, long long b) {
  return a == 0 ? b : gcd(b % a, a);
}

// return (n / m)
long long jacobi(long long n, long long m) {
  if (n == 1)
    return 1;
  if (n % 2 == 0)
    return ((m*m-1)/8 % 2 == 0 ? 1 : -1) * jacobi(n/2, m);
  return ((n-1)/2*(m-1)/2 % 2 == 0 ? 1 : -1) * jacobi(m % n, n);
}

// n が素数か判定
bool solovay_strassen(long long n) {
  if (n % 2 == 0) return n == 2;  
  random_device seed_gen;
  default_random_engine engine(seed_gen());
  uniform_int_distribution<> dist(1, n - 1);
  for (int i = 0; i < 20; ++i) {
    long long a = dist(engine);
    if (gcd(a , n) > 1) return false;
    if ((jacobi(a, n) + n)% n != pow(a, (n - 1) / 2, n))
      return false;
  }
  return true;
}

参考文献

Solovay, Robert, and Volker Strassen. “A fast Monte-Carlo test for primality.” SIAM journal on Computing 6.1 (1977): 84-85.

平方剰余の相互法則がここまでダイレクトにアルゴリズムの設計に根差したものは初めて見ました。賢い。