カメヲラボ

主にプログラミングとお勉強全般について書いてます

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

数学

はじめに

素数の問題は小学生くらいの知識であっても取り組むことが出来る。しかし奥深さは世界の天才数学者をも悩ませる。
ドイツが生んだ大数学者ガウスは言った。

「素数は数学のょぅι゛ょであるﾊｧﾊｧ」

うそ。言ってない。それはさておき、コンピュータと素数に魅了されてしまった私(と素敵な素数野郎たち)が、2年以上の歳月をかけて悩み・戦い・築き上げたシュゴー（゜Д゜）なテクニックを惜しげもなく晒し、素数初心者・プログラミング初級者にもできるだけわかりやすく解説した「素数マスターへの道」を、突然ではありますが連載したいと思う次第でございます。ょぅι゛ょは出てきません。

前置きはいいから

さて、何をするのかというと、以下のサイト
http://www.spoj.pl/problems/PRIC/
にある問題に対して正しい答えを出力する、超高速なプログラムを作ろう！というのが主旨です。

「Prime Checker」という問題で、内容はいたってシンプル。

1から始めて、そこに1234567890をどんどん加えていく。ただし、2^31(2の31乗=2147483648)を超えると、2^31で割って「余り」を使う。そうやって計算した数が、素数ならば「1」素数で無いならば「0」と出力しなさい。

いくつか具体的に書くと、

1
1234567891
321652133(1234567891+1234567890は2^31を超えるので余りを取る)
1556220023
643304265(1556220023+1234567890は2^31を超えるので余りを取る)
1877872155
964956397(1877872155+1234567890は2^31を超えるので余りを取る)
52040639(964956397+1234567890は2^31を超えるので余りを取る)
1286608529
373692771(128660852+1234567890は2^31を超えるので余りを取る)
...

のような数が、素数なのか素数でないのかを判定していきます。どこまで調べるかというと、1から数えて33,333,333個(3千万個くらい)について、です。

とりあえずやってみるのをオススメします

これから数学的・プログラミング的なテクニックを紹介・解説していきますが、本当に楽しむには自分でやってみるのが一番ですので、最初に述べておきます。もちろん、読むだけで十分楽しめるようには工夫しようと思いますが。

素数かどうかを調べる基本的なお話

その数が素数であるかどうかを確かめるには、実際に割り算してみるのが一番単純な方法です。たとえば7が素数かどうかは、2,3,4,5,6と割り算してみて「確かに割り切れる数が無い！」と確認すればよいわけです。ただし、本当に間の数全部で試す必要は無くて、もし「2で割り切れなければ」4とか6みたいな「2の倍数ではすべて割り切れない」ということがはっきりしています。ですので、試しに割ってみるといっても、すべての数ではなく、「間にある素数」で割るだけで確かめられます。もっと言えば、大きすぎる数は商が0になってしまうので、その数の平方根(2乗してその数になる数)まで試せば十分です。

もうちょっと大きな数で言えば、67が素数かどうかを調べるには、8*8=64ですから、8までの素数、2,3,5,7で割れるかを調べれば良いことになります。

いきなり壁が

問題の紹介でいくつかの数を書きましたが、1の次がいきなり1234567891になっています。12億ちょっとです。10桁あるので、普通の電卓に収まりきらない、かなり大きな数です(;´Д｀)

幸い私達が使っているパソコンには簡単に計算する機能があるので、たとえばWindowsなら電卓のアプリが付いています。「1234567891」と打ち込んでから「x^y」というボタンを押し、そのあと「0.5」入れて「=」ボタンを押します。すると、
35136.41830067487114928194813943
のように表示されました。小数点以下は切り捨てておいて、35136までの素数で割り切れるかどうかを調べれば良いことになります。しかしこれは1234567891の例ですので、最大でいくつになるのかを調べておきます。扱うのは2^31を超えない数ですので、一番大きい数は2^31より1小さい数の2147483647です。これの平方根を計算すると、約46340になりました。つまり、事前に46340までの素数リストを作る必要があります。

エラトステネスの篩い(ふるい)

素数のリストを作る方法として最も有名な方法です。具体的なことは、ウィキペディア
http://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9
で非常に分かりやすく説明されています。要は、2から始めて倍数をどんどん消していく方法で、プログラムを書くのも簡単です。コードはそこいらじゅうに溢れていますが、やっぱり自分で書かないとダメですね。

#include<stdio.h>
#include<stdlib.h>
#include<string.h>

#define N 46340

int genPrime(char* check, int*primeList )
{
  int nPrime=0;
  int i,j;
  
  for(i=2;i<N;++i)
  {
    if(check[i]==0)
    {
      for(j=i+i;j<N;j+=i)check[j]=1;
      primeList[nPrime]=i;
      nPrime+=1;
    }
  }

  return nPrime;
}

int main()
{
  char check[N]={0}; /* 0なら素数1なら素数ではない */
  int primeList[N]={0};
  int nPrime; /* 見つかった素数の数 N=46340ならnPrimeは4792になるはず */
  int i;
  
  nPrime = genPrime(check, primeList);
  
  for(i=0;i<nPrime;++i)
  {
    printf("%d ",primeList[i]);
  }
  printf("\n");
  printf("%d prime num.\n",nPrime);

  return 0;
}

46340までの素数は4792個あることがわかったので、一つの数が素数かどうかを判定するには最大4792回試し割をすれば良いわけです。回数はちょっと多い気もしますが、最近のコンピュータなら一瞬で計算してくれます。

本当に一瞬か？

何にしても、素数のリストが出来たのでとりあえずプログラムを書いてみます。

/* 0.1.c */
/* 割り算で素数判定 */
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#define bool unsigned char
#define true 1
#define false 0

#define N 46340
#define D 1234567890U
#define R 2147483648U

/* M個の数について素数判定をする */
#define M 33333333

int genPrime(char* check, unsigned int*primeList )
{
  int nPrime=0;
  int i,j;
  
  for(i=2;i<N;++i)
  {
    if(check[i]==0)
    {
      for(j=i+i;j<N;j+=i)check[j]=1;
      primeList[nPrime]=i;
      nPrime+=1;
    }
  }

  return nPrime;
}

bool isPrime(unsigned int*pList, int nPrime, unsigned int n)
{
  bool is_prime = (n!=1); /* 1の時だけfalseにしておく */
  int i;

  for( i=0; pList[i]<n && i<nPrime; ++i)
  {
    if(n%pList[i]==0)
    {
      is_prime=false;
      break;
    }
  }
  
  return is_prime;
}

int main()
{
  char check[N]={0}; /* 0なら素数1なら素数ではない */
  unsigned int primeList[N]={0};
  int nPrime; /* 見つかった素数の数 N=46340ならnPrimeは4792になるはず */
  int i;
  unsigned int n=1;
  
  /* 最初に素数のリストを作っておく */
  nPrime = genPrime(check, primeList);
  
  for(i=0;i<M;++i)
  {
    if(isPrime(primeList, nPrime, n)) printf("1");
    else printf("0");
    n=(n+D)%R;
  }
  printf("\n");
  
  return 0;
}

これをgccでコンパイル(-O3オプションで最適化)して実行してみると、AMD Sempron 3000+(1.79GHz)でおよそ540秒(9分)かかりました。言い忘れましたが、SPOJでは制限時間が25秒で、サーバマシンの性能はそれほど高くないので、一応解けるものの制限時間内には全然収まりません。SPOJのサーバにpostしてみたところ、75万個くらいまでは判定できましたが、全部を判定するには少なくとも45倍くらいスピードアップしなければなりません。
さてさて、どないしましょ。

続く