次にアルゴリズムであるが、いろいろいじってみた結果こうした。

1. まず p < 256 な小さな素数でエラトステネスのふるいにかけ
2. 次にMiller-Rabin素数判定法を適用する

これは「個々の64bit整数が素数かどうか」を判定するのには(素数表を引くことを除けば)ベストに近い方法なのだが、「任意の範囲の整数の素数を全て上げよ」という問題には実は適していない。「任意の範囲」を全てメモリー上に置けるのであれば、そこに直接エラトステネスのふるいをかけた、残った数をあげた方が速いのだ。

まずは32bitの素数表をあらかじめ作っておく。

mksieve32.c

/*
 * $Id: mksieve32.c,v 0.1 2010/07/27 15:21:58 dankogai Exp dankogai $
 */
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>

#include <stdint.h>
typedef uint32_t U32;
typedef uint64_t U64;

#include "bitmap.c"

int main(void) {
    bitmap *b = bitmap_new((U64)INT_MAX+1, NULL);
    bitmap_fill(b, 1);
    bitmap_set( b, 0, 0); /* 1 is not prime */
    U64 p, i;
    for (p = 3; p < USHRT_MAX ;){
        printf("sieving %llu\r", p);
        fflush(stdout);
        for ( i = p + p ; i <=  UINT_MAX ; i += p ) {
            if ((i & 1) == 0) continue; /* skip even */
            bitmap_set( b, i>>1, 0 );
        }
        for(p += 2; !bitmap_get(b, p>>1); p += 2); /* seek next prime */
    }
    char filename[]  = "sieve32.bm";
    printf("saving %s\n", filename);
    bitmap_save(b, filename);
    return 0;
}

これ自体、(もうすぐ一世代遅れになる?)Core i7なiMacだと30秒ほどしかかからない。なお、bitmap.cを拡張して、bitmap_fill()とbitmap_save()を追加した。これらが何をするかは自明であろう。ちなみにbitmap_load()がないのは、bitmap_new()にその機能がすでに存在するから。

primes.c

あとはこの32bit素数表を利用すれば原理的には64bit整数の素数を全て数え上げることができるわけだ。ただし今度は全整数にふるいをかけるのではなく、あくまで範囲指定した範囲のみにかける。その点に留意すれば、mksieve32.cとやってることは大して変わらない。

/*
 * $Id: primes.c,v 0.2 2010/07/27 15:35:59 dankogai Exp dankogai $
 */
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <errno.h>
#include <math.h>

#include <stdint.h>
typedef uint32_t U32;
typedef uint64_t U64;

#include "bitmap.c"

#define SIEVEFILE "sieve32.bm"

int main(int argc, char **argv) {
    bitmap *sieve = bitmap_new(0, SIEVEFILE);
    if (!sieve){
        perror(SIEVEFILE);
        exit(errno);
    }
    U64 start = argc > 1 ? atoll(argv[1]) : 0;
    U64 size  = argc > 2 ? atoll(argv[2]) : 1000*1000*1000;
    bitmap *b = bitmap_new(size >> 1, NULL);
    bitmap_fill(b, 1);
    if (start == 0) bitmap_set(b, 0, 0);
    U64 pmax = (U64)sqrtl(start+size);
#ifdef VERBOSE
    printf("pmax = %llu\n", pmax);
#endif
    U64 p, i;
    for (p = 3; p < pmax;){
#ifdef VERBOSE
        printf("sieving %llu\r", p);
        fflush(stdout);
#endif
        for ( i = p + p - (start % p) ; i <= size ; i += p ) {
            if ((i & 1) == 0) continue; /* skip even */
            bitmap_set( b, i>>1, 0 );
        }
        for(p += 2; !bitmap_get(sieve, p>>1); p += 2); /* seek next prime */
    }

    char filename[256];
    snprintf(filename, 256, "%llu~%llu.bm", start, start+size);
#ifdef VERBOSE
    printf("saving %s\n", filename);
#endif
    bitmap_save(b, filename);
    return 0;
}

これを使うと、10兆以下であれば10秒弱で範囲10億の検索が終わる。999,999,999,000,000,000~1,000,000,000,000,000,000、すなわち10京近傍の10億ですら15秒ほどであった。10兆であれば特に並列化しなくとも1日程度で終わるということになる。必要なディスク容量は61Mの一万倍なので610Gといったところか。

量子コンピュータ
竹内 繁樹

記者が最初にプログラムを書いたのは、電気店の店頭にあった8ビット機だが、当時のマイコンではとてもこんなところまでは来れなかった。できるだけ長生きして、またチャレンジしてみたいと思う。

いずれは「サマーウォーズ」にもちらっと登場したShorのアルゴリズムを実装した量子コンピュータが登場して、RSAに使われているような大きな数の素数判定も一瞬で終わるようになる日が来るかもしれない。しかしその日が来ても、エラトステネスのふるいはさまざまな形で使われ続けるだろう。もしかして、人類滅亡後も。

Dan the Prime Counter

bitmap.c

/*
 * $Id: bitmap.c,v 0.2 2010/07/27 15:27:02 dankogai Exp dankogai $
 */
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <fcntl.h>
#include <sys/stat.h>
#include <sys/mman.h>

typedef struct {
    int    fd;
    size_t size;
    char  *map;
} bitmap;

bitmap *bitmap_free(bitmap *b){
    if (b){
        if (b->map) {
            if (b->fd){
                munmap(b->map, b->size);
                close(b->fd);
            }else{
                free(b->map);
            }
        }
        free(b);
    }
    return (bitmap *)NULL;
}

bitmap *bitmap_new(size_t size, const char *filename){
    bitmap *b = (bitmap *)malloc(sizeof(bitmap));
    struct stat st;
    if (!b) return (bitmap *)NULL;
    if (filename){
        b->fd = open(filename, 
                     size ? O_RDWR|O_CREAT : O_RDONLY, 
                     size ? (mode_t)0644   : (mode_t)0444);
        if (b->fd == -1) {
            perror(filename);
            return bitmap_free(b);
        }
        if (size){
            if (ftruncate(b->fd, size >> 3) == -1){
                perror(filename);
                return bitmap_free(b);
            }
            b->map = (char *)mmap(0, size >> 3, PROT_READ|PROT_WRITE,
                                  MAP_SHARED, b->fd, 0);

        }
        else{
            fstat(b->fd, &st);
            size = st.st_size << 3;
            b->map = (char *)mmap(0, size >> 3, PROT_READ, 
                                  MAP_PRIVATE, b->fd, 0);
        }
        if (b->map == MAP_FAILED) return bitmap_free(b);
    }
    else{
        if (!size) return (bitmap *)NULL;
        b->map = (char *)malloc(size >> 3);
        if (!b->map) return bitmap_free(b);
        b->fd = 0;
    }
    b->size = size;
    return b;
}

int bitmap_save(bitmap *b, const char *filename){
    int ok = 0, fd = open(filename, O_RDWR|O_CREAT|O_TRUNC, (mode_t)0644);
    if (fd != -1){
        if (write(fd, b->map, (b->size >> 3)) != -1) ok = 1;
    }
    if (!ok) perror(filename);
    close(fd);
    return ok;
}

void bitmap_fill(bitmap *b, int val){
    size_t i;
    for (i = 0; i < (b->size >> 3); i++) b->map[i] = val ? 0xff : 0;
}

static const int bits[] = { 1, 2, 4, 8, 16, 32, 64, 128 };

inline int bitmap_set(bitmap *b, size_t where, int val){
    if (val) b->map[where >> 3] |=  bits[where & 7];
    else b->map[where >> 3] &= ~bits[where & 7];
    return val;
}

inline int bitmap_get(bitmap *b, size_t where){
    return !!(b->map[where >> 3] & bits[where & 7]);
}

#ifdef TEST
#include <errno.h>
int main (int argc, char **argv){
    bitmap *b = bitmap_new(4096, (argc > 1 ? argv[1] : NULL));
    if (!b) return errno;
    bitmap_set(b, 0, 1);
    printf("%d\n", bitmap_get(b, 0));
    bitmap_free(b);
    return 0;
}
#endif