/* * 以下の記事のプログラム http://www.gizmodo.jp/2016/03/prime-number-not-so-random.html > 素数の出方はランダムではなかった。1億個調べて浮かんだ奇妙な数 元の論文 http://arxiv.org/abs/1603.03720 *著作権 Copyright(c) 2016 松元隆二。 *ライセンス GPL。 * 使い方 このファイルをコンパイルしてください。64bitOSでのみ確認。 gcc count-prime.c -O3 -Wall **素数データのダウンロード 以下のA or B, C どれかの素数データをダウンロードしてください。 A. ftp://ftp.ai.kyutech.ac.jp/Numbers/prime/2-32.diff.txt.gz 最大素数2^32までの素数データ(diff形式) 512MByte B. ftp://ftp.ai.kyutech.ac.jp/Numbers/prime/10-12.diff.txt.gz 最大素数10^12までの素数データ(diff形式) 30GByte C. ftp://ftp.ai.kyutech.ac.jp/Numbers/prime/10-11.diff.txt.gz 最大素数10^11までの素数データ(diff形式) 10GByte 巨大なファイルですのでmd5でのパリティチェックをおすすめします。 ftp://ftp.ai.kyutech.ac.jp/Numbers/prime/md5sum.txt パリティチェック このデータ(diff形式)は以下のようなデータです。 ---- 2 1 2 2 (以下続く) ---- 最初の素数は2で、次の素数を求めるには次の数字を足してください。 ---- 2 -> 素数2 1 -> 2+1 -> 素数3 2 -> 3+2 -> 素数5 2 -> 5+2 -> 素数7 ---- このデータ形式にしているのは、ファイルのサイズ削減のためです。 **計算実行。 引数無しで実行すると以下のように表示されます。 > $ ./a.out > usage: ./a.out: 進数 [素数の数] > 進数は3以上。 > 標準入力に差分形式の素数表を入力 引数で進数を指定してください。 オプション引数で素数の数を指定してください。未指定の場合は 標準入力に入ったデータすべてを対象とします。 論文1ページ目の3進数、π(x0) = 10^6 の計算をする場合は $ gzip -dc hoge-diff.txt.gz | ./a.out 3 100000 としてください。 -- 以下の出力が表示されると思います。 --- 標準入力データ読み飛ばし素数=2 標準入力データ読み飛ばし素数。2ペア。=(2, 3) 標準入力データ読み飛ばし素数。2ペア。=(3, 5) 集計結果 論文に合わせたPI(x0)=1000000 本当のPI(x0)=1000001 最後に対象とした素数=15485867 PI(x0;3,0)=1 PI(x0;3,1)=499829 PI(x0;3,2)=500170 本当のPI(x0)=1000003 最後に対象とした素数=15485927 PI(x0;3,(1,1))=215873 PI(x0;3,(1,2))=283957 PI(x0;3,(2,1))=283957 PI(x0;3,(2,2))=216213 --- 最初に素数をいくつか読み飛ばすため、論文でπ(x0) = 10^6と書いてあって も実際に計算対象となる素数はそれ以上になります。それを本当のPI(x0)と して表示します。 なお、第2引数を省略した場合は、読み込んだ素数の数を表示するだけで 本当のPIとか論文に合わせたPIとか表示しません。 ------- $ gzip -dc 2-32.diff.txt.gz | ./a.out 3 標準入力データ読み飛ばし素数=2 標準入力データ読み飛ばし素数。2ペア。=(2, 3) 標準入力データ読み飛ばし素数。2ペア。=(3, 5) 集計結果 読み込んだ素数の数 PI(x0)=203280221 最後に対象とした素数=4294967291 PI(x0;3,0)=1 PI(x0;3,1)=101638091 PI(x0;3,2)=101642128 最後に対象とした素数=4294967291 PI(x0;3,(1,1))=45563873 PI(x0;3,(1,2))=56074218 PI(x0;3,(2,1))=56074218 PI(x0;3,(2,2))=45567909 ------ 論文にある10進数, 10^8までの結果は以下のようにします。 -------- $ gzip -dc 2-32.diff.txt.gz | ./a.out 10 100000000 標準入力データ読み飛ばし素数=2 標準入力データ読み飛ばし素数=3 標準入力データ読み飛ばし素数。2ペア。=(2, 3) 標準入力データ読み飛ばし素数=5 標準入力データ読み飛ばし素数。2ペア。=(3, 5) 標準入力データ読み飛ばし素数=7 標準入力データ読み飛ばし素数。2ペア。=(5, 7) 標準入力データ読み飛ばし素数。2ペア。=(7, 11) 集計結果 論文に合わせたPI(x0)=100000000 本当のPI(x0)=100000004 最後に対象とした素数=2038074793 PI(x0;10,1)=24999437 PI(x0;10,3)=25000135 PI(x0;10,7)=25000401 PI(x0;10,9)=25000027 本当のPI(x0)=100000005 最後に対象とした素数=2038074803 PI(x0;10,(1,1))=4623042 PI(x0;10,(1,3))=7429438 PI(x0;10,(1,7))=7504612 PI(x0;10,(1,9))=5442345 PI(x0;10,(3,1))=6010982 PI(x0;10,(3,3))=4442562 PI(x0;10,(3,7))=7043695 PI(x0;10,(3,9))=7502896 PI(x0;10,(7,1))=6373981 PI(x0;10,(7,3))=6755195 PI(x0;10,(7,7))=4439355 PI(x0;10,(7,9))=7431870 PI(x0;10,(9,1))=7991431 PI(x0;10,(9,3))=6372941 PI(x0;10,(9,7))=6012739 PI(x0;10,(9,9))=4622916 -------- 以上です。 */ #include #include #include #include #include #include int main (int argc, char **argv) { int debug_input = 1; //入力データの簡単なチェック if (argc < 2 || argc > 4) { printf ("usage: %s: 進数 [素数の数]\n", argv[0]); printf (" 進数は3以上。\n"); printf (" 標準入力に差分形式の素数表を入力\n"); exit (1); } int base = strtol (argv[1], NULL, 10); if (base < 3) { printf ("ERROR: 異常な進数=%d\n", base); exit (1); } uint64_t max_line = UINT64_MAX; if (argc == 3) { max_line = strtoull (argv[2], NULL, 10); if (max_line < 1) { printf ("ERROR: 異常な素数の数=%" PRIu64 "\n", max_line); exit (1); } } uint64_t count1[base]; uint64_t count2[base][base]; memset (count1, 0, base * sizeof (uint64_t)); memset (count2, 0, base * base * sizeof (uint64_t)); uint64_t line_count = 0; uint64_t num = 0, pre_num = 0; uint64_t line_count_sub1 = 0; uint64_t line_count_sub2 = 0; uint64_t line_count_sub1_save = 0; uint64_t line_count_sub2_save = 0; uint64_t count1_last = 0; uint64_t count2_last = 0; while (!feof (stdin)) { if (line_count_sub2 == max_line && line_count_sub1 == max_line) break; const int BUFFSIZE = 1024; char buf[BUFFSIZE + 1]; buf[0] = '\0'; fgets (buf, BUFFSIZE, stdin); if (buf[0] == '\0') break; if (buf[0] == '\n') continue; if (buf[0] == '#') continue; uint64_t diff = strtoul (buf, NULL, 10); if (diff == 0) { printf ("ERROR: 異常な入力=\"%s\"\n", buf); exit (1); } line_count++; pre_num = num; num += diff; if (line_count <= 4 && debug_input) { /*簡単な入力テスト。 最初の入力データは 2 -> 素数2 (line1) 1 -> 2+1 -> 3 (line2) 2 -> 3+2 -> 5 (line3) 2 -> 5+2 -> 7 (line4) のはず。 */ int flag = 0; switch (line_count) { case 1: flag = (num == 2); break; case 2: flag = (num == 3); break; case 3: flag = (num == 5); break; case 4: flag = (num == 7); break; } if (flag == 0) { printf ("ERROR: 標準入力データのフォーマットエラー=\"%s\"\n", buf); exit (1); } } // 元論文に合わせた。 if (num < base) { printf ("標準入力データ読み飛ばし素数=%" PRIu64 "\n", num); if (pre_num != 0) { printf ("標準入力データ読み飛ばし素数。2ペア。=(%" PRIu64 ", %" PRIu64 ")\n", pre_num, num); } continue; } if (line_count_sub1 < max_line) { line_count_sub1++; count1[num % base]++; line_count_sub1_save = line_count; count1_last = num; } // 元論文に合わせた。 if (pre_num <= base) { printf ("標準入力データ読み飛ばし素数。2ペア。=(%" PRIu64 ", %" PRIu64 ")\n", pre_num, num); continue; } if (line_count_sub2 < max_line) { line_count_sub2++; count2[pre_num % base][num % base]++; line_count_sub2_save = line_count; count2_last = num; } } if (max_line != UINT64_MAX) { if (line_count_sub2 != max_line || line_count_sub1 != max_line) { printf ("ERROR: 標準入力データの素数の数が足りません\n"); exit (1); } } //出力 printf ("\n集計結果\n"); if (max_line == UINT64_MAX) { printf ("読み込んだ素数の数 PI(x0)=%" PRIu64 "\n\n", line_count); } else { printf ("論文に合わせたPI(x0)=%" PRIu64 "\n\n", max_line); } if (max_line != UINT64_MAX) { printf ("本当のPI(x0)=%" PRIu64 "\n", line_count_sub1_save); } printf ("最後に対象とした素数=%" PRIu64 "\n\n", count1_last); int i; for (i = 0; i < base; i++) { uint64_t n = count1[i]; if (n == 0) { continue; } printf ("PI(x0;%d,%d)=%" PRIu64 "\n", base, i, n); } printf ("\n"); if (max_line != UINT64_MAX) { printf ("本当のPI(x0)=%" PRIu64 "\n", line_count_sub2_save); } printf ("最後に対象とした素数=%" PRIu64 "\n\n", count2_last); for (i = 0; i < base; i++) { int j; for (j = 0; j < base; j++) { uint64_t n = count2[i][j]; if (n == 0) { continue; } printf ("PI(x0;%d,(%d,%d))=%" PRIu64 "\n", base, i, j, n); } } return 0; }