10兆以下の素数の一覧を求める

いささか旧聞に属するのですが、IT Proに10兆までの素数のリストを作ってみませんか？という記事があります。2010年5月の記事なので2年近く前ですね。この記事に触発されて、色々な人がこの命題に挑戦していたようですが、恥ずかしながら私はそれを知りませんでした。今さらですが私も挑戦してみたいと思います。

単純に10兆以下の素数一覧を求めるというだけであれば、素数生成プログラムを作成し、ひたすら実行し続ければいつかは終了することでしょう。ただし、それではこのお題を達成したことになりません。ただ単にプログラムを作って実行するだけではなく、その見積もりも行うことが含まれています。つまり、かなりビジネスよりの命題ということですね。とはいえ、10兆以下の素数を単純な方法で求めようとすると計算時間がとても現実的ではありませんから、このような見積もりは事前に必要になるところです。

なお、IT Proの記事は2ページ構成になっていて後半は会員登録(無料)をしないと読めません。私は会員登録をしていないので、後半は読んでいません。記者がどのように見積もり、どのような結果を得たかは知らずに以下の文章とプログラムは書かれています。この場合、結果を読まずにやった方が楽しめていいかもしれませんけれどね(^^)。

では具体的に見積もりの検討に入りましょう。まず想定できるのは、10兆という数値は32ビット整数の範囲には収まらないということです。これは単純に64ビット整数を用いるようにしましょう。記者は32ビット整数に収まる範囲では32ビット整数を用い、それ以上では64ビット整数を用いるという方法をとったそうですが、私は面倒なので全部64ビット整数で処理することにします。このような場合分けを行わなければならないほど、今時のPCのメモリ搭載量は少なくないですし、場合分けを行う開発・実行のコストも無視できないと考えます。

次に検討が必要なのは、結果データの格納方法です。素数一覧によると、10億以下の素数一覧データはzip圧縮して125MBになります。これは展開すると490MBにもなります。10兆以下で単純に1万倍になるとして5TBが必要になってしまいます。PC用のHDDも2012年時点で単体で4TBのものが販売されていますし、RAIDを組むことで5TBのファイルを格納することも出来ますが、それでは確かに芸がない。とりあえず素数をそのまま出力するのではなく、直前に出力した素数との差分で出力するようにプログラムを改良したところ、10億以下で137MBとなりました。これで想定ファイルサイズは1.5TBとなりかなり現実的になってきました。更に、記者が行ったようにデータを圧縮することで10億以下が36MBになりました。10兆以下だと350GB程度の見積もりになるでしょう。通常ならこれで問題ないサイズなのですが、残念ながら私の使っているPCのHDDに空きが無いため、これでも収まりません。ということで、ここは「実現は可能だ」という結論を得たことで満足し、実行時にはデータを保存しないことにします。

次に実行時間の見積もりです。ちなみに単純な方法をとった場合にどれくらい掛かるでしょうか。素数で書いている通り、1億以下で5分、10億以下で2時間かかりました。探索範囲が10倍に対して計算時間が24倍と仮定したとして、10兆以下なら663,552時間掛かることになります。これは約76年ですからとても現実的とは言えませんね。記者も採用しているエラトステネスのふるいを用いることで、1億以下が4秒、10億以下が44秒です。探索範囲が10倍に対して計算時間が13倍と仮定して10兆以下の計算時間は約7日と現実的な時間に収まりそうです。

実行に使用したプログラムは以下のものです。基本的に素数で作成したプログラムを64bit化しただけです。

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

#define ARRAY_SIZE_MAX   (100000000)
#define MAX_N (10000000000000)

char    array[ARRAY_SIZE_MAX];
long long primelistcount;
long long *primelist;

int initprimelist( void );
long long isqrt( long long x );

int main( int ac, char *av[] )
{
    long long base, arraysize, n, i, j, prev;

    /*  序数側の素数一覧を生成  */
    if( !initprimelist() )
        return( 0 );

    prev = 0;
    /*  ARRAY_SIZE_MAX分ごとの整数区間をふるいにかける   */
    for( base = 2; base < MAX_N; base += ARRAY_SIZE_MAX ) {
        /*  整数区間配列の大きさを決める    */
        if( MAX_N - base + 1 > ARRAY_SIZE_MAX )
            arraysize = ARRAY_SIZE_MAX;
        else
            arraysize = MAX_N - base + 1;

        /*  配列を初期化する    */
        for( i = 0; i < arraysize; i++ )
            array[i] = 1;

        /*  配列をふるいにかける    */
        for( i = 0; i < primelistcount; i++ ) {
            n = primelist[i];
            if( base <= n )
                j = n;
            else if( base % n == 0 )
                j = base;
            else
                j = ( base / n + 1 ) * n;
            j = j < n * n ? n * n : j;
            for( ; j < base + arraysize; j += n )
                array[j - base] = 0;
        }

        /*  ふるいで残った数は素数である    */
        for( n = 0; n < arraysize; n++ ) {
            if( array[n] == 1 ) {
                printf( "%lld\n", base + n - prev );
                prev = base + n;
            }
        }
    }

    return( 0 );
}

/*  序数側の素数一覧を生成  */
int initprimelist( void )
{
    long long n, i;
    long long sqrt_llong_max;

    /* sqrt(MAX_N)を求める */
    sqrt_llong_max = isqrt( MAX_N );

    /*  配列を初期化する    */
    for( i = 0; i <= sqrt_llong_max; i++ )
        array[i] = 1;

    /*  配列をふるいにかける    */
    for( n = 2; n <= sqrt_llong_max; n++ ) {
        if( array[n] == 1 ) {
            for( i = n * n; i <= sqrt_llong_max; i+=n ) {
                array[i] = 0;
            }
        }
    }

    /* 素数一覧配列にコピー */
    primelistcount = 0;
    for( n = 2; n <= sqrt_llong_max; n++ ) {
        if( array[n] == 1 )
            primelistcount++;
    }
    primelist = calloc( primelistcount, sizeof(long long) );
    if( primelist == NULL )
        return( 0 );
    for( n = 2, i = 0; n <= sqrt_llong_max; n++ ) {
        if( array[n] == 1 )
            primelist[i++] = n;
    }

    return( 1 );
}

long long isqrt( long long x )
{
    long long s, t;

    if( x == 0 ) return 0;
    s = 1;
    t = x;
    while( s < t ) {
        s <<= 1;
        t >>= 1;
    }
    do {
        t = s;
        s = ( x / s + s ) >> 1;
    } while( s < t );

    return( t );
}

見積もりの検証も兼ねて、10以下から10倍ずつ実行してみました。1兆以下なら1日以内に終了する見積もりですので、それほど時間をとられるわけではなく問題ありません。

1兆以下の素数を求めるのに約18時間、圧縮後のデータサイズ約38GBですので、ここまで大体見積もり通りとなっています。ということで、いよいよ10兆以下の素数を求めますが、上にも書きましたとおり結果を格納するだけのHDD空き容量がありませんのでwcコマンドでデータサイズだけを計測することにします。その結果は以下の通りとなりました。

実行時間

392,993.64秒

素数の数

346,065,536,839個

データサイズ

972,251,362,068バイト

実行時間が予想より少々短いのが気になりますが、その他の結果は大体予想通りですので、おそらく問題ないでしょう。HDDの空き容量が確保できたら、改めて検証してみたいと思います。

(2015/12/3追記)

今頃になってですが、放置していたことに気が付いたので検証してみました。マシンが変わっているので参考ですが、実行時間は582,396.92秒、圧縮後データサイズは292,202,856,977バイトでした。素数の数はもちろん一致しています。

(2016/9/2追記)

10兆以上の一覧を出力するには手元のストレージが不足するのですが、個数を調べるだけなら出来るかな。ということで100兆以下の素数の個数を求めてみましたところ、約3か月かかって3,204,941,750,802個という結果が出ました。

(2020/4/25追記)

ソースコードを少し整理しました。