10兆以下の素数の一覧

2012/4/3作成

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 MAX_ARRAY   (100000000)
#define MAX_N (10000000000000)
#define START_R (0x40000000)

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

int initprimelist( void );

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

    if( !initprimelist() )
        return( 0 );

    prev = 0;
    /*  MAX_ARRAY分ごとの整数区間をふるいにかける   */
    for( base = 2; base < MAX_N; base += MAX_ARRAY ) {
        /*  整数区間配列の大きさを決める    */
        if( MAX_N - base + 1 > MAX_ARRAY )
            arraysize = MAX_ARRAY;
        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, r;
    long long sqrt_llong_max;

    /* sqrt(MAX_N)を求める */
    r = START_R;
    sqrt_llong_max = 0;
    while( r ) {
        if( MAX_N >= ( sqrt_llong_max + r ) * ( sqrt_llong_max + r ) )
            sqrt_llong_max += r;
        r /= 2;
    }

    /*  配列を初期化する    */
    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 );
}

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

範囲 実行時間(秒) 素数の数 データサイズ(バイト) 圧縮後データサイズ(バイト)
10以下 0.00 4 8 33
100以下 0.00 25 50 49
1000以下 0.00 168 368 138
1万以下 0.00 1229 2878 746
10万以下 0.00 9592 23622 5601
100万以下 0.03 78498 200066 46917
1000万以下 0.39 664579 1734013 423374
1億以下 4.28 5761455 15307793 3894399
10億以下 51.43 50847534 137066738 36044465
100億以下 583.58 455052511 1241646236 335816423
1000億以下 6178.82 4118054813 11356614904 3147200553
1兆以下 66116.03 37607912018 104717120836 38042643653

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

実行時間(秒)
392993.64
素数の数
346065536839
データサイズ(バイト)
972251362068

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

(2015/12/3追記)

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

(2016/9/2追記)

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


あおやぎのさいと2.0初歩の整数論プログラミング