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

フェルマーテストの実装

まずは、基数を色々変えられる一般的なものを書いてみます。

/* a^n mod n を求める */
unsigned int mpow(unsigned int a, unsigned int n)
{
  unsigned int buf[32]={0}; /* a^(2^k) mod n を保持 */
  unsigned int ans=1, r=n;
  int k=0;
  
  buf[0]=a; /* a^(2^0)==a */
  
  for(k=1; (1<<k)<n; ++k)
  {
    buf[k] = ( (unsigned long long)buf[k-1] * buf[k-1] ) % n;
    assert(k<32);
  }
  
  for(; k>=0; --k)
  {
    if( (1<<k)<=r )
    {
      ans = ( (unsigned long long)ans * buf[k] ) %n;
      r-=(1<<k);
    }
  }
  assert(r==0);
  
  return ans;
}

32ビット整数同士の乗算部分は倍精度にして、桁溢れを防いでおきましょう。

フェルマーテストで再度挑戦！

基数を2として、フェルマーテストによる素数判定プログラムを実行してみます。先ほどのmpowという関数を以下の様に使用しました。

#define D 1234567890U
#define R 2147483648U
#define M 33333333

int main()
{
  int i;
  unsigned int n=1;
  
  for(i=0;i<M;++i)
  {
    if( mpow(2,n)==2 )printf("1");
    else printf("0");
    n=(n+D)%R;
  }
  printf("\n");
  
  return 0;
}

実行時間は初回と同じ環境(http://d.hatena.ne.jp/Ozy/20100225#p1)でおよそ70秒でした。約7.7倍!かなり速くなりましたね。しかしフェルマーテストは、前回で述べたように不完全です。33333333個のうち、242個を誤判定してしまいました。そこで、基数を3としたテストも重ねて行ってみます。

    if( mpow(2,n)==2 )printf("1");

の部分を

    if( mpow(2,n)==2 && mpow(3,n)==3 )printf("1");

と書き換えて実行すると、8秒程度遅くなりましたが、誤判定の数は53個まで減りました。3以外の基数で試してみると、

 5 -> 50
 7 -> 44
11 -> 54

の様な結果になったので、組み合わせとしては2と7が良さそうです。

続く