10兆までの素数リスト
プログラミングコンテストチャレンジブックを読んでいたら、素数リストの作り方が出てきた。
そういえば、１０兆までの素数リストを作ってみる、なんて記事があったな。
これだ
http://itpro.nikkeibp.co.jp/article/Watcher/20100519/348242/?ST=develop
読み返しているうちに、ついつい挑戦してみる気になってしまった。
色々試して意外だったのが、10億までの素数をリストアップした場合でも「試し割り」より「エラトステネスの篩」の方が100倍以上も高速だったことだ。リストが大きくなればその差はますます開いていく。
メモリ上の広い範囲をまんべんなくアクセスするエラトステネスの篩は、実際には効率が悪いのではないかと思ったりしたのだが、見当が外れてしまった。
考察してみるに
- 試し割りは意外と無駄が多い。まず素数の密度だがn/log(n)に従うらしい。100億あたりだとおよそ5%だった。大雑把にlog(n)の項を無視して「ほぼnに比例する」と考える。そして試し割りだが、これまた大雑把にいえば、3について無駄になる確率は2/3, 5については(2/3)*(4/5)、7は(2/3)*(4/5)*(6/7)、11は(2/3)*(4/5)*(6/7)*(10/11)だ。√nまで続くから、nまで探索したときの試し割りの回数はO(n*√n)で抑えられると言えそうだ。一方、エラトステネスの篩は3についてはn/3回、5についてはn/5回といった処理回数が必要だから、全体ではn*{(1/3)+(1/5)+(1/7)+...}回になる。これは少なくともO(n*log(n))で抑えられる。正確なことはよく分からないが要するに計算オーダーが違うようだ。
- 64bit整数の除算が遅い。エラトステネスの篩では基本的に除算がない。（ただし探索区間を分割した際に手間を惜しんで使用している）
- 一度に篩にかける範囲を500KByte程度に制限することで、CPUのキャッシュの恩恵を受けることができた
といった理由が考えられるだろうか。
「試し割り一部→ミラー・ラビン→試し割り完全」も試してみたが100億程度の数では効果は微かで、エラトステネスの篩には全く敵わなかった。ミラー・ラビンは単一の数を素数判定するには有効だが、リストを作るのには向かないということだろうか。
（ミラー・ラビンの実装は http://d.hatena.ne.jp/kazu-yamamoto/20100131/1264921432 を参考にした）
そんなわけで「エラトステネスの篩」を基本に素数リスト生成プログラムを作成した。以下概要。
- まずふるい落とすべき素数のリストを「試し割り」で作成する。10兆まで欲しいならば、10兆^(1/2)、つまり316万程度までのリストを用意すればよい（合成数であれば、必ず√n以下の素因数をもつ）
- CPUのキャッシュ(実験環境ではそれぞれのコアにL2キャッシュ=512KB)を考慮して、一度に篩にかける範囲を50万以下とし、これを繰り返した。
- 実験環境は４コアのCPU（AMD Phenom II X4 945, 3.0GHz）だったので、マルチスレッドによる分散処理を行った。スレッドを３つ作成し、それぞれのスレッドは探索ブロックを３つおきに処理する。ブロックが小さいため、これらをバラバラのファイルに出力するのはよくない。そこで同期処理し各スレッドが担当したブロックを順番に出力するようにした。幸い、各スレッドの処理時間はほぼ等しいため、同期してもCPU使用率は300%に近い値になった。高速化の効果は2.1倍くらいになった。スレッドを４つではなく３つにしたのは、出力結果をパイプで受け取り圧縮するプロセスが１つあることに配慮したものだ。４つスレッドをたてた場合は、プロセスの切り替えが多発するためか、かえって性能が悪化してしまった。
- 2*3*5*7*11*13*17=510,510周期の「17までふるい落し済みリスト」を使い回す方針も試したが100億まででさえ10%程度の高速化の効果しか得られず、それ以上になればさらに効果が薄れること、実装が多少複雑になることから最終的には採用しなかった。
さて問題は出力だが、あいにく余り容量の大きいディスクを持っていない。
差分をVariableByte形式にして、さらにgzipで圧縮してみたものの、それでも10兆までのリストを保持するには約230GBが必要という見積もりになった。
3以降の素数と素数の差分は必ず偶数で2以上なので、「2で割って1を引く」という処理を加えてみたが、当然、雀の涙の効果しかなかった。
どうしたらいい…買うか？いやいや、そんな。
そうだ…保存しなければいいんだ！
目的は処理時間、およびファイルサイズを知ることだ。その結果230GBのファイルを保存できたとしても、すぐに消去することになるだろう。それならば最初から保存しなければよい。
time ./prime | gzip -c > primes.dat
と
time ./prime | gzip -c | wc -c
を比較して、処理時間が殆ど同じことを確かめた。ディスクIOはボトルネックではない。
そんな訳で、めでたく10兆までの素数リストを作成するのに必要な時間とデータサイズが見積もれた。
| どこまでか | 処理時間 | サイズ（gzip圧縮） |
|---|---|---|
| 10億まで | 3秒45 | 約30.8MB |
| 100億まで | 28秒89 | 約287MB |
| 1000億まで | 4分19秒 | 約2.68GB |
| 1兆まで | 51分48秒 | 約25.1GB |
| 10兆まで | 13時間27分 | 約236GB |
プログラム
#include <stdio.h> #include <stdlib.h> #include <inttypes.h> #include <cmath> #include <algorithm> #include <pthread.h> uint64_t g_check_max = 1000ULL*1000*1000*1000*10; uint64_t g_search_width = 1000ULL*500; uint64_t* g_primes; int g_thread_count = 3; pthread_mutex_t g_mutex = PTHREAD_MUTEX_INITIALIZER; volatile uint64_t g_seq_number = 0; volatile uint64_t g_last_out = 3; struct ThreadInfo{ pthread_t thread; int thread_id; uint64_t from; }; void out_prime(uint64_t v){ if(v <= g_last_out) return; uint64_t d = (v - g_last_out) / 2 - 1; g_last_out = v; while(d > 0x80){ fputc(0x80 | (d & 0x7f), stdout); d >>= 7; } putc(d, stdout); } bool isPrime(uint64_t v){ if(v <= 1) return false; if(v == 2) return true; if((v&1) == 0) return false; for(uint64_t i = 3; i * i <= v; i+=2){ if(v % i == 0) return false; } return true; } void sift_prime(uint64_t block_number, char* ar, uint64_t from, uint64_t to){ uint64_t size = to - from; for(uint64_t i = 1; g_primes[i]*g_primes[i] < to; i++){ uint64_t p = g_primes[i]; for(uint64_t j = from + p - 1 - ((from + p - 1) % p); j < to; j += p){ ar[j-from] = 0; } } while(true){ pthread_mutex_lock(&g_mutex); if(g_seq_number == block_number) break; pthread_mutex_unlock(&g_mutex); usleep(100); } for(uint64_t i = 1 - (from & 1); i < size; i+=2){ if(ar[i]) out_prime(from+i); } g_seq_number++; pthread_mutex_unlock(&g_mutex); } void search_range(ThreadInfo* info){ uint64_t width = g_search_width; char* ar = new char[width]; uint64_t block_number = info->thread_id; for(uint64_t from = info->from + width * info->thread_id; from < g_check_max; from += width * g_thread_count){ uint64_t to = std::min(g_check_max, from + width); memset(ar, 1, width); sift_prime(block_number, ar, from, to); if(block_number / g_thread_count % 100 == 0){ fprintf(stderr, "thread:%d, %lld : %lld\n", info->thread_id, from, to); } block_number += g_thread_count; } } int main(){ uint64_t total = 0; uint64_t check_max_r = (uint64_t)sqrt(g_check_max) + 1; g_primes = new uint64_t[check_max_r]; g_primes[0] = 2; uint64_t primes_i = 1; uint64_t i; for(i = 3; true; i += 2){ if(isPrime(i)){ out_prime(i); g_primes[primes_i] = i; primes_i++; if(i >= check_max_r){ i += 2; break; } } } ThreadInfo* threads = new ThreadInfo[g_thread_count]; for(int th = 0; th < g_thread_count; th++){ threads[th].thread_id = th; threads[th].from = i; pthread_create(&(threads[th].thread), NULL, (void*(*)(void*))search_range, &threads[th]); } for(int th = 0; th < g_thread_count; th++){ pthread_join(threads[th].thread, NULL); fprintf(stderr, "thread end:%d\n", th); } return 0; }
付録＊
■シングルスレッドとマルチスレッドの比較
100億まで
コマンドは time ./prime | gzip | wc -c
| シングルスレッド | 3スレッド | 4スレッド |
|---|---|---|
| 65秒51 | 28秒89 | 38秒01 |
■一度に処理する探索範囲の大きさの比較
100億まで
| スレッド | 範囲 | 処理時間 |
|---|---|---|
| シングルスレッド | 5万個 | 60秒77 |
| シングルスレッド | 50万個 | 65秒51 |
| シングルスレッド | 500万個 | 129秒74 |
| 3スレッド | 5万個 | 161秒67 |
| 3スレッド | 50万個 | 28秒89 |
| 3スレッド | 500万個 | 87秒16 |
■結果を保存するかどうかの比較
100億まで、3スレッド、VariableByte
| コマンド | 処理時間 |
|---|---|
| time ./prime | gzip -c > primes.dat | 29秒37 |
| time ./prime | gzip -c > /dev/null | 28秒56 |
| time ./prime | gzip -c | wc -c | 29秒18 |
| time ./prime > /dev/null | 23秒64 |
■出力データの大きさの比較
100億まで、3スレッド、ファイルに保存する
| 方式 | 処理時間 | サイズ(Byte) |
|---|---|---|
| テキスト | 93秒56 | 4,948,214,537 |
| テキスト+gzip | 4分45秒 | 1,160,824,945 |
| 差分テキスト | 60秒02 | 1,241,646,236 |
| 差分テキスト+gzip | 3分28秒 | 353,280,158 |
| 差分VariableByte | 24秒05 | 455,615,134 |
| 差分VariableByte+gzip | 33秒12 | 287,681,964 |
| 差分VariableByte(x/2-1) | 24秒49 | 455,052,953 |
| 差分VariableByte(x/2-1)+gzip | 29秒21 | 286,813,303 |
テキスト+gzipの場合、gzipがボトルネックになったようだ。
double数値のテキストシリアライズとバイナリシリアライズの速度比較
uchiumikさんが、CRFの拡張のようなものを作成しているのだが、パラメーターが多くなりすぎて、モデルデータの読み込みに時間がかかっているという。何百万、何千万個のdouble値をファイルから読み込む処理をしているという。
そのせいで実験がしんどくては、結局手を抜いても特をしないので、そこは是非改善すべきだと思うのだが、説得力を出すために実験してみた。（本当は自分がやってみたかっただけだったりする）
本当は共有メモリをオススメしたいのだが、シリアライズをテキストからバイナリに変更するだけでもそれなりに効果があるはずだ。
次のコードで実験
#include <stdio.h> #include <list> void write_bin1(FILE* fp, double* data, int size){ fwrite(data, 1, sizeof(double)*size, fp); } void read_bin1(FILE* fp, double* data, int size){ fread(data, 1, sizeof(double)*size, fp); } void write_text(FILE* fp, double* data, int size){ for(int i = 0; i < size; i++){ fprintf(fp, "%e\n", data[i]); } } void write_text_prec(FILE* fp, double* data, int size){ for(int i = 0; i < size; i++){ fprintf(fp, "%.15e\n", data[i]); } } void read_text(FILE* fp, double* data, int size){ const int BUF_SIZE = 1000; char buf[BUF_SIZE]; for(int i = 0; i < size; i++){ fgets(buf, BUF_SIZE, fp); data[i] = atof(buf); } } void read_text_scan(FILE* fp, double* data, int size){ for(int i = 0; i < size; i++){ fscanf(fp, "%lf", data + i); } } int main(){ const int test_size = 1000000; FILE* fp; double* data = (double*)malloc(sizeof(double)*test_size); for(int i = 0; i < test_size; i++){ data[i] = rand() * (1.0 / RAND_MAX); } //* fp = fopen("test.dat", "wb"); write_text_prec(fp, data, test_size); fclose(fp); /*/ fp = fopen("test.dat", "rb"); read_text(fp, data, test_size); fclose(fp); //*/ free(data); return 0; }
(コンパイラはgcc version 4.0.1、コンパイル方法は g++ -O2 -DNDEBUG binary_io.cpp)
結果
| 方法 | write | read |
|---|---|---|
| テキスト16桁 | 2.282sec | 1.026sec(2.189sec) |
| テキスト7桁 | 0.739sec | 0.248sec(0.441sec) |
| バイナリ | 0.098sec | 0.046sec |
double精度を保持するためには、仮数部は十進で16桁程度必要とのことで、16桁で実験した。
7.826369259425611e-06 1.315377881431662e-01 7.556053221950332e-01 4.586501319234493e-01
のような形で保存される。精度を指定しなかった場合は7桁で保存されるようだ。
もちろん、バイナリの場合はこのような心配は不要。
実験の結果、バイナリの方が約20倍以上高速ということになった。（データの初期化のみの処理時間は0.028secなので、その分を差し引いた方がより正確になるか）
もちろん、デメリットにも注意する必要はあると思う。
データファイルがエンディアン依存、doubleの形式に依存してしてしまうことと、ファイルを人間が読むのがつらいというのが主な点だろうか。
表のカッコの中は、fscanf()を使った場合の結果で、テキストファイルを扱う場合でも、fgets() + atof()を使ったほうが良いことが分かる。
ファイルがキャッシュされているかで結果が変わるかもしれないが、それは同等の条件ということにしておこう。
もう一つ、この部分にこだわっても、ボトルネックじゃなかったりするとあまり嬉しくないので、その点を確認しておく。
ファイルから読み取ったあと、std::list
void write_bin1(FILE* fp, double* data, int size){ fwrite(data, 1, sizeof(double)*size, fp); } void read_bin1(FILE* fp, double* data, int size){ fread(data, 1, sizeof(double)*size, fp); } void write_bin2(FILE* fp, const std::list<double> data){ for(std::list<double>::const_iterator it = data.begin(); it != data.end(); it++){ fwrite(&(*it), 1, sizeof(double), fp); } } void read_bin2(FILE* fp, std::list<double>* data, int size){ data->resize(0); for(int i = 0; i < size; i++){ double val; fread(&val, 1, sizeof(double), fp); data->push_back(val); } } void write_bin3(FILE* fp, const std::list<double> data){ double* buf = (double*)malloc(sizeof(double) * data.size()); int i = 0; for(std::list<double>::const_iterator it = data.begin(); it != data.end(); i++, it++){ buf[i] = *it; } write_bin1(fp, buf, data.size()); free(buf); } void read_bin3(FILE* fp, std::list<double>* data, int size){ double* buf = (double*)malloc(sizeof(double) * size); read_bin1(fp, buf, size); data->resize(0); for(int i = 0; i < size; i++){ data->push_back(buf[i]); } free(buf); }
結果は
| 方法 | write | read |
|---|---|---|
| bin1 | 0.098sec | 0.046sec |
| bin2 | 0.551sec | 0.492sec |
| bin3 | 0.530sec | 0.474sec |
bin3の方法は、freadやfwriteの呼び出し回数を減らす工夫をしたのだが、bin2と殆ど変わらなかった。
もちろん、このlistのコストは、テキストファイルの場合も追加で発生する。
そこから予測すると、listに入れる場合はバイナリ化しても2〜4倍程度しか高速化は期待できない。
蛇足だが、perlとpythonでもバイナリを読み込んでみた
open(IN, '<', 'test.dat'); binmode(IN); read(IN, $indata, 1000000*8); @data = unpack('d*', $bindata);
0.291sec (perl 5.8.8)
import struct fp = open('test.dat', 'rb') bindata = fp.read(1000000*8) data = struct.unpack('1000000d', bindata)
0.202sec (python 2.5.1)
結構はやい
でも実験結果を考えると、使い方によっては、それほどdouble値のテキストがボトルネックにはならないこともありそうだ(このままいくと1000万個でも23秒程度だ)。バイナリにすることで不便になる点もあるので、慎重に検討したい。
一閃突きの確率
いまさらだけどドラクエ９の話題。
「一閃突き」とは、1/2の確率でミス、1/2の確率でかいしんの一撃が出るヤリのスキルで、メタル系のスライムにも有効なので重宝するワザだ。
オノのスキルの「魔人切り」も同じ効果だが、こちらはアニメーションが若干長いため、一閃突きの方が好まれていると思われる。
さて、このワザだが、「ミスが出るときは連続してミスが出る」という噂がある。「一閃突き」は本当にミスが連続しやすいのだろうか？
ミスの確率が1/2だとすると、理論上、5回連続でミスする確率はたったの1/32になる。
とはいえ、それを目撃するために平均5*32=160回試さないといけないわけじゃない。
5回以上連続で試したとすると、1〜5回目と2〜6回目と3〜7回目・・・
と連続5回になるチャンスがじゃんじゃん増える。
そんなわけで
( 1/32 ) + ( ( n - 5 ) * ( 1/64 ) ) = 1
n=67
で,67回で5連続以上の回数の期待値が1回になる。
（6連続以上のとき5連続以上を1回とカウントするために、2回目以降から検出するときは直前がヒットだった場合のみカウントするので1/32ではなく1/64で計算）
67回の一閃突き、つまり13.4回の戦闘で、5回連続ミスを1回くらい目撃する
132回の一閃突き、つまり26.4回の戦闘で、6回連続ミスを1回くらい目撃する
261回の一閃突き、つまり52.2回の戦闘で、7回連続ミスを1回くらい目撃する
・・・
予想に反して、連続ミスは起こりやすいのかもしれない。
さて、前置きはともかく、実験してみた。
- ゴールデンスライムなどの敵を相手にひたすら一閃突きをする
- 戦闘が終了しても、連続してカウントを続ける。MPの回復などは適宜行う
- 「かいしんの一撃」と表示されたら、かわされても成功とカウントする
- 試行回数659回（飽きたところで終了）
この実験は、敵の行動など乱数の消費がどうという考察はなくアバウトなものである点、ご了承頂きたい。一応、8割方パーティー4人が先制攻撃したと思うので、ミスの連続に偏りがあれば、傾向が出てくるはずだと思われる。
以下ではヒットを1、ミスを0で表記する
- ミスとヒットの回数
0:276
1:283
- 2連続のパターンの回数
00:127
01:149
10:148
11:134
- 3連続のパターンの回数
000:57
001:70
010:83
011:66
100:69
101:79
110:65
111:68
- 4連続のパターンの回数
0000:25
0001:32
0010:39
0011:31
0100:42
0101:41
0110:35
0111:30
1000:32
1001:37
1010:44
1011:35
1100:27
1101:38
1110:30
1111:38
- 5連続のパターンの回数
00000:13
00001:12
00010:12
00011:20
00100:21
00101:18
00110:15
00111:15
01000:19
01001:23
01010:21
01011:20
01100:17
01101:18
01110:16
01111:14
10000:12
10001:20
10010:26
10011:11
10100:21
10101:23
10110:20
10111:15
11000:13
11001:14
11010:23
11011:15
11100:10
11101:20
11110:14
11111:24
以上、実験結果からは、著しい偏りは見られなかったと言えるでしょうか。
下記はpythonの実験コードです。
hit = "0010100001110100010000111000110101000110101010010101100100001001001011101011110100110010010100111011101000011001111110010101100000011110111111100100000011101111010101111101000101011111110000001101100010001101110110100000100010110110101111111110011100000111101100011100010011011001110101010000111101101010101001101000000001011000001101101001001111000110101111010001001001011101010011100010101010111010111010001100101100010010110001101101100110010010010001100100101010100100011101001001011010010110010110100101100101111111111101110100110110010001000111111100011" for w in (1,2,3,4,5): freq = {} for i in xrange(len(hit)-w+1): key = hit[i:i+w] if key not in freq: freq[key] = 0 freq[key] += 1 for (key,val) in sorted(freq.items(), lambda x, y: cmp(x[0], y[0])): print "%s:%d" % (key,val)
■
はてなダイアリー、開設カンタンでいいなあ。
貧乏性なので、カスタマイズしてタイトル部分の高さを節約してみた。
ブログのタイトルは「雷雲の変奏曲」
今日は雷鳴がすごかったので・・・。
プログラミングの話題を中心に書き進めていきたいと思います。
最近コツコツ取り組んでいるのがコンピューターオセロ。
世界で最初に8x8オセロを解けたらいいな。
まあムリだろうけど。
6x6は解かれているらしい。
同一局面数（鏡像・回転も統合）を数えてみた
| 着手 | 盤面数 | 終局盤面数 | 増率 |
|---|---|---|---|
| - | 1 | 0 | - |
| 1手目 | 1 | 0 | 1.00 |
| 2手目 | 3 | 0 | 3.00 |
| 3手目 | 14 | 0 | 4.67 |
| 4手目 | 60 | 0 | 4.29 |
| 5手目 | 322 | 0 | 5.37 |
| 6手目 | 1,773 | 0 | 5.51 |
| 7手目 | 10,649 | 0 | 6.01 |
| 8手目 | 67,245 | 0 | 6.31 |
| 9手目 | 434,029 | 36 | 6.45 |
| 10手目 | 2,958,586 | 28 | 6.82 |
| 11手目 | 19,786,652 | 925 | 6.69 |
| 12手目 | 137,642,678 | 1,321 | 6.96 |
| 13手目 | 912,004,217 | 26,938 | 6.63 |
| 14手目 | 6,158,977,944 | 46,489 | 6.75 |
| 15手目 | 38,755,255,016 | 578,729 | 6.29 |
| 16手目 | 244,248,103,534 | ? | 6.41 |
今のところ数えることが出来たのはここまで。(2010/07/17 更新) (2013/07/27 表だけ更新14手→16手)
普通に考えると、同一局面を纏めるためには、纏められる可能性のある局面をすべてメモリにキャッシュしておかなければならない。
1盤面あたり16byte程度必要なので、6,158,977,944通りでは約96GBのメモリが必要になる。
そんなにメモリを積んだマシンは持っていないので、メモリを節約する工夫をしている。
盤面に対してハッシュ値を計算したのである。例えばハッシュ値を0〜40とすることで、ハッシュ値が0のものだけ保存して同一盤面かチェック、残りは捨てる。次はハッシュ値が1のもの、2のもの・・・と41回繰り返すことで、処理時間は約41倍になるが、メモリの使用量は約1/41にすることができた。こうして4GBのメモリで足りるようになった。
実際の処理時間はおよそ36時間だった。
オセロってあまりやらないけど、10手目で約300万ものバリエーションがあるものなのかと感心してしまう。
この実験では、石が四辺に到達することで、着手の分岐数が少なくなり、同一局面数の増加が次第に穏やかになるのではないか？という仮説を検証したかった。
しかし、14手目まででは、まだちょっとよく分からない。
15手目まで求めるには、このやり方では約7*7=49倍の時間がかかる。あまり現実的ではない。
