新年早々頭の体操、inspired by Team-lablog » 1000万個目の素数を超高速に出力せよ。
素数とは、言うまでもなく、中学校の数学でちょっと頭を痛める「1以外の、自分自身と1以外に約数がない正の自然数」だ。
どうやったら求められるだろうか？
コンピュータは基本的にバカだが力持ちだ。
我々がやりたくないような、「1からその数までを割る数にして、順番に割り算をして割り切れるか確認する」などというアホな方法でも求めてくれる。
しかしこれでは時間が掛かりすぎる。
ちょっと頭をひねれば、
- 2以外の偶数は素数ではないので、はなから計算しない
- その数以外に割り切れる数があれば、そこで計算をやめる
とすれば、計算量が少なくて済む=早く計算が終わることが分かる。
しかし、この方法では素数かどうか確認する数の大きさに比例して計算量が増えていく。
// This is C# source code
ulong count = 1;
ulong seed;
ulong target;
bool isSosu;
target = 3;
while()
{
isSosu = true;
for(seed = 2; seed < target; count++)
{
if(count % seed == 0)
{
isSosu == false;
break;
}
}
if (isSosu) count++;
if (count >= 1000000) break; // stop the loop when 1Mth sosu is found;
target += 2;
}
100までとか、1000まで、位なら良いが、2^32 = 4,294,967,296まで、とかやると、非常に単純な(偶数を省かない、途中で割り切れても計算を止めない場合)計算量はざっと(N-1) + (N-2) + ... 1だ、死ぬる。
さて、素数を求めるアルゴリズムで一番よく知られているのが、「エラトステネスの篩」である。
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
これは、N以下の素数を求めるのには
(1)sqrt(N)までの素数を求め、
(2)それらの素数の倍数で無いものが素数
という解法だ。
まず(1)について、例えば100までの数の場合には、sqrt(100) = 10以下の素数、すなわち2、3、5、7の倍数でなければ、素数であるということである。
ここからちょっと発展させると、2^(2^n)以下の数は2^(2^n-1)より小さい素数の倍数でなければ、素数ということになる。
2^(2^5) = 2^32 = 4,294,967,296までの数は、2^(2^4) = 2^(16) = 65,536までの素数の倍数でなければ素数、ということになる。
単純計算なら4,294,967,296/1から4,294,967,296/4,294,967,296まで計算(4,294,967,296回)しなくてはらないところを、4,294,967,296/1から4,294,967,296/65,536まで計算する(65,536回)だけでよい、ということになれば相当減ったのが分かるだろう。
しかも既知の素数でしか計算しないので、上記の場合でも65,536回よりはぐっと少なくなる(おおよそ1/10)はずである。
さらに上記の方法と同様に偶数を計算しないようにすれば、更に速くなるはずだ。
(2)は電子計算機特有の問題の解になる。
PCも含め、電子計算機は一般的に乗算よりも除算が遅い。
このため「割り切れるか」を見るよりは、「素数の倍数か」どうかを確認する方が、速度が速くなる。
なので、以下のように篩を実装した。
(補記 on 20150115、C#からCに変更)
(補記 on 20150119、ソースを整理)
pIsSosuという配列を作り、奇数番目のものはPRIMARY_NUMBER = 1で、偶数番目のものはNOT_PRIMARY_NUMBER = 0で初期化しておく。
// This is C code
#include "stdafx.h"
#include "stdlib.h"
#include "time.h"#define PRIMARY_NUMBER 1
#define NOT_PRIMARY_NUMBER 0
#define MAX_NUMBER 268435456 // 16384^2 = (2^14)^2
#define SEED_MAX 16384
#define COUNT_SOSU_MAX 10000000int _tmain(int argc, _TCHAR* argv[])
{
clock_t clock_start = clock();
// print time
printf("Start the 10Mth primary number calculation\n");
//
unsigned long seed, count;
unsigned long countSosu;
char* pIsSosu = (char*)malloc(sizeof(char)* MAX_NUMBER);
// print time
printf("The filter array allocated: %dmsec \n", clock() - clock_start);
// initialize array, all even numbers are flagged as non primary
for (count = 0; count < MAX_NUMBER - 1; count += 2)
{
*(pIsSosu + count) = PRIMARY_NUMBER;
*(pIsSosu + count + 1) = NOT_PRIMARY_NUMBER;
}
*pIsSosu = NOT_PRIMARY_NUMBER; // for 1
*(pIsSosu + 1) = PRIMARY_NUMBER; // for 2
// print time
printf("The filter array initialized: %dmsec \n", clock() - clock_start);
//
countSosu = 1; // 2 is counted
for (seed = 3; seed < SEED_MAX;)
{
// count up Sosu
countSosu++;
// flag the multiple of seed as "not primary number
for (count = seed * seed - 1; count < MAX_NUMBER; count += seed)
{
*(pIsSosu + count) = NOT_PRIMARY_NUMBER;
}
// find next primary number from pIsSosu array
seed += 2; // even numbers are not primary, so skip any time
while (*(pIsSosu + seed - 1) == NOT_PRIMARY_NUMBER) seed += 2;
}
// print time
printf("The primary number array obtained: %dmsec \n", clock() - clock_start);
//
for (count = SEED_MAX; count < MAX_NUMBER; count++)
{
if (*(pIsSosu + count) == PRIMARY_NUMBER) countSosu++;
if (countSosu >= COUNT_SOSU_MAX) break;
}
// print time and result
printf("The target obtained: %dmsec, %d\n", clock() - clock_start, count + 1);
printf("Press Return key to finish.\n");
getchar();
// clean-up heap
free(pIsSosu);
return 0;
}
筆者のPCのメモリーの関係で、pIsSosuが最大2^32 / 8 = 536,870,912までしか確保することができなかった。
これでもおそらく1000万番目の素数までは十分数えられるだろう。
pIsSosuのindexは1ずらしてあり、1が0、2が1...という風になっている。
最初に1に相当する*(pIsSosu)はNOT_PRIMARY_NUMBERに、2に相当する*(pIsSosu + 1)はPRIMARY_NUMBERにしておく。
以下forループでseedを3からスタートし、seedの倍数に相当する*(pIsSosu + index)をNOT_PRIMARY_NUMBERに設定していく。
そして次の素数を探し、同じ作業を続けていく。
さてさて結果はというと、あっという間に100万番目の素数、15,485,863を見つけることができた。
ただ何か早くないんだよね(10秒強)、要検討。
(追記)
お、早くなった。
100万番目までなら、最大数が2^24 = 16,777,216、確認する素数の範囲がsqrt(2^24) = 2^12 = 4,096まででも十分なので、これらを変えれば1秒以内になりました、めでたし。
(追記 on 20150115)
ソースをCに書き換えて、VS2013でrelease buidで試したところ、
- 100万番目: 0.43秒
- 1000万番目: 7.22秒
もう少し何とかなるのかな。
(追記 on 20150119)
ソースを上記のものに変えて、Sandy Bridge Core i5-2500で回したところ、トータル3.18秒、ふるいの部分だけなら2.8秒、もうここが限界か？
(追記 on 20150119その2)
ソースを一部変更、素数の倍数をpIsSosuにフラグしていく際に、開始点は2 * seed - 1ではなくて、seed * seed - 1でよい (seed^2未満の数は既にそれより小さいseedのループでフラグされている)。
これでトータル3.13秒。
