素数マスターへの道(4)

すげーぜ！ミラー・ラビン

フェルマーテストによって、高確率で素数判定が出来るようになりましたが、もうちょっと確率の高い判定法があります。
ミラー・ラビン（とかラビン・ミラー）判定法と呼ばれる方法を使うと、同じくらいの時間で、誤判定をさらに減らすことが出来ます。
nが素数かどうかを判定するために、フェルマーテストでは

a^n mod n ≡ a

が成り立つかどうかを調べました。

ミラー・ラビンテストでは、n(奇数に限る)を使って、以下のように行います。
(1)n-1は偶数なので、

n-1 = 2^k * q

の形に変形できます。

(2)(1)で求めたqを使って、

a^q mod n ≡ 1

かどうかを調べます。これが成り立つ場合、nが素数である可能性が非常に高いと言えます。

(3)(2)が成り立たない場合は、kより小さいi{0,1,2,...,k-1}について、

(a^q)^(2^i) mod n ≡ n-1

が成り立つかどうかを調べます。上記の式が成り立つiが1つでもあれば、nは素数である可能性が非常に高いと言えます。

(4)(2)と(3)の判定をすべて通り抜けた場合は、素数ではありません。

フェルマーテストと比べると若干複雑ではありますが、計算方法自体は大きく変わるものではないので、実装が難しいというほどではありません。

bool mrtest(unsigned int a, unsigned int n)
{
  unsigned int buf[32]={0}; /* a^(2^k) mod n を保持 */
  unsigned int ans=1, r;
  unsigned int q=n-1;
  int l=0, k=0, j=0, i=0;
  
  /* n-1 = 2^j * q の形にする */
  j=0;
  for(; q!=0 && q%2==0; )
  {
    ++j;
    q/=2;
  }
  
  /* この辺はmpow関数と同じ */
  buf[0]=a; /* a^(2^0)==a */
  for(k=1; (1<<k)<q; ++k)
  {
    buf[k] = ( (unsigned long long)buf[k-1] * buf[k-1] ) % n;
    assert(k<32);
  }
  
  /*  まずは a^q mod n の判定 */
  r=q;
  l=k;
  for(; l>=0; --l)
  {
    if( (1<<l)<=r )
    {
      ans = ( (unsigned long long)ans * buf[l] ) %n;
      r-=(1<<l);
    }
  }
  assert(r==0);
  if(ans==1||ans==n-1)return true;
  
  /* (a^q)^(2^i) mod n の判定 */
  for(i=1; i<j; ++i)
  {
    ans = ( (unsigned long long)ans * ans ) %n;
    if(ans==n-1)return true;
  }
  
  return false;
}

これを用いて実行すると、基数を2とした場合59個の誤判定がありました。フェルマーテストの場合と同じように、2と7の組み合わせでやると誤判定は2個だけ！
2個のうち1つはn=1の場合なので、簡単に処理できますし、実質1つと考えて良いでしょう。この1つは、例えば2,5,7のような基数の組み合わせで3回テストすれば正しく判定できるのですが、そうすると結構時間がかかってしまうので、予め誤判定する値をコードの中に埋め込んでしまうのが良さそうです。59個しか誤判定がないのであれば、テストする基数も2だけにしておいても良さそうですね。

続く