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

シングルスレッドでは５万個の場合が最もよいが、3スレッドでは５０万個が最もよいという結果になった。一度の処理が小さすぎると、スレッドの同期処理が負担になるからだと考えられる。

■結果を保存するかどうかの比較
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

テキスト形式の場合、HDDがボトルネックになったようだ。
テキスト+gzipの場合、gzipがボトルネックになったようだ。