世間は夏休みだそうだし、連日の猛暑で体調も底だし、というわけで私も素数を数えてみた。

10兆までの素数のリストを作ってみませんか？ - 記者の眼：ITpro
もしあなたがプログラマだったら、プログラムを書いて10兆までの素数のリストを作ってみてほしい。情報システムの開発に携わる人であれば、10兆までの素数のリストを出力するシステムの見積もりを考えてみてほしい。費用はどれくらいかかるか、納期はどれくらいか、あなたはどんな答を出すだろうか。仕様書はうまく書けるだろうか。
asin:4873113806
プライムナンバーズ
David Wells / 伊知地宏監訳 / さかいなおみ訳
[原著:Prime Numbers: The Most Mysterious Figures In Math]

といっても原田記者と同じように書いても芸がないので、同記事につっこみつつ、私はどうかいたのかを解説していくことにする。

10兆までの素数のリストを作ってみませんか？ - 記者の眼：ITpro
次の面白さは、大きな出力結果に対処することだ。図2と図3では出力結果を文字列型（String）の変数に連結しているが、この方法では10億まで探索すると現実的な時間では終わらない。10億までの素数は（私が計算したところでは）5084万7534個あり、それをシンプルなテキストファイルにすると527Mバイトの大きさになる。この大きな出力をどう扱うか、プログラマの基本的な技能が試される。

すぐに思いつくのはビットマップで結果を保存することだろう。これだとUINT_MAX、232までの素数表も4GBの1/8、512MiBに収まることになる。さらに付け加えれば、素数は2を除けば全て奇数なので、奇素数のみを保存するようにすればさらに半分、256MiBに収まる。高々10億であれば高々61MB程度である。

というわけでbitmap.cというものをこさえた。ビットマップを簡単に扱うためのライブラリである。bitmap_new()でファイル名を指定すると、そこにmmap()するので、書き出しと圧縮は不要である。その代わり結果を読み出すには、bits2txt.cを使う。

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

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

このMiller-Rabin素数判定法、一般的には確率的素数判定法に用いられるが、「小さな」素数であれば決定的素数判定法として使えるのだ。

Miller-Rabin primality test - Wikipedia, the free encyclopedia
When the number n to be tested is small, trying all a < 2(ln n)2 is not necessary, as much smaller sets of potential witnesses are known to suffice. For example, Pomerance, Selfridge and Wagstaff[6] and Jaeschke[7] have verified that
  • if n < 1,373,653, it is enough to test a = 2 and 3;
  • if n < 9,080,191, it is enough to test a = 31 and 73;
  • if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
  • if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
  • if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
  • if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
cover cover
素数入門/数論入門
芹沢正三

要するに、 341兆5500億7172万8320 以下であれば、底は 2, 3, 5, 7, 11, 13, 17 までチェックするだけでよいということになる。

このMiller-Rabinをはじめ、数論の世界では冪剰余という演算がよく出てくる。冪乗の段階で超巨大な数字が出てくるので大変かと思いきや、剰余を取りながら冪乗することにより意外に効率よく計算できるのが面白い。

ところがその効率のよいアルゴリズムでも、オーバーフローが生じてしまう。これをどうするかだが、powmod.cでは二通りの解決法を用意した。GCCの128bit演算サポートを利用する方法と、GMPを使う方法だ。GCC 4.4なUbuntu 10.04 LTSでは前者が、GCC 4.2なSnow Leopardでは後者の方が高速だった。

というわけで、primes.cである。Core i7 な iMac で、10億の範囲を検索するのに1プロセス300秒前後、5分といったところである。x86_64なUbuntuとOS Xで動作検証してある。

素数を数えたいかたはご自由にご利用くださいませ。

Dan the Prime Counter

bitmap.c

/*
 * $Id: bitmap.c,v 0.1 2010/07/26 07:00:51 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;
}

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

bits2txt.c

Compile:
gcc -Wall -W -O6 -o bits2txt bits2txt.c 
#include <stdio.h>
#include <errno.h>
#include "bitmap.c"

int main(int argc, char **argv){
    if (argc < 2){
	fprintf(stderr, "%s filename_of_bitmap\n", argv[0]);
	exit(-1);
    }
    bitmap *b = bitmap_new(0, argv[1]);
    unsigned long long offset = atoll(argv[1]);
    if (!b) {
	perror(argv[1]);
	exit(errno);
    }
    size_t i;
    for (i = 0; i < b->size; i++){
	if (bitmap_get(b, i)) printf("%llu\n", offset + (i<<1) + 1);
    }
    return 0;
}

powmod.c

/*
 * $Id: powmod.c,v 0.1 2010/07/26 07:00:51 dankogai Exp dankogai $
 */

#ifdef USE128BIT

typedef unsigned int U128 __attribute__((mode(TI)));
U64 powmod(U64 base, U64 power, U64 mod){
    U64 result = 1;
    while (power > 0){
        if (power & 1) 
            result = result <= UINT_MAX && base <= UINT_MAX 
                ? (result * base) % mod
                : ((U128)result * (U128)base) % mod;
        base = base <= UINT_MAX 
            ? (base * base) % mod
            : ((U128)base * (U128)base) %mod;
        power >>= 1;
    }
    return result;
}

#else

#include <gmp.h>

inline U64 powmod_gmp(U64 base, U64 power, U64 mod){
    mpz_t b, p, m, r;
    mpz_init(b);
    mpz_init(m);
    mpz_init(p);
    mpz_init(r);
    mpz_set_ui(b, base);
    mpz_set_ui(p, power);
    mpz_set_ui(m, mod);
    mpz_powm(r, b, p, m);
    return mpz_get_ui(r);
}


inline U64 powmod_64(U64 base, U64 power, U64 mod){
    U64 result = 1;
    while (power > 0){
        if (power & 1) result = result * base % mod;
	base = (base * base) % mod;
	power >>= 1;
    }
    return result;
}

inline U64 powmod(U64 base, U64 power, U64 mod){
    return base  >= UINT_MAX ? powmod_gmp(base, power, mod)
	:  power >= UINT_MAX ? powmod_gmp(base, power, mod)
	:  mod   >= UINT_MAX ? powmod_gmp(base, power, mod)
	:  powmod_64(base, power, mod);
}

#endif

primes.c

Compile:
gcc -Wall -W -O6 -I/opt/local/include -L/opt/local/lib -lgmp -o primes primes.c
or
gcc -Wall -W -O6 -DUSE128BIT -o primes primes.c
Usage:
primes [start=0] [size=1,000,000,000]
/*
 * $Id: primes.c,v 0.1 2010/07/26 07:00:51 dankogai Exp dankogai $
 */
#include 
#include 
#include 

typedef int Bool;
#define FALSE 0
#define TRUE  1

typedef unsigned long long U64;

#include "powmod.c"


/* primes up to 256 */
static U64 seed_primes[]  = {
    /*0,1,2,3, 4, 5, 6, 7 */
    2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
    101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
    191,193,197,199,211,223,227,229,233,239,241,251, 0
};

/*
 * http://en.wikipedia.org/wiki/Miller–Rabin_primality_test
 * works for n < 341,550,071,728,321
 */
#define MRMAX 7 /* till 17 */

static inline Bool millar_rabin(U64 a, U64 s, U64 d, U64 n){
    U64 r;
    U64 x = powmod(a, d, n);
    if (x == 1) return TRUE;
    for (r = 0; r < s; r++){
	if (x == n - 1) return TRUE;
	x = (x * x) % n;
    }
    return FALSE;
}


Bool is_prime(U64 n){
    U64 d = n - 1, s = 0, i;
    if (n <= 1) return FALSE;
    if ((n & 1) == 0) return n == 2;
    for (i = 1; seed_primes[i]; i++) {
        if (n == seed_primes[i]) return TRUE;
        if (n % seed_primes[i] == 0) return FALSE;
    }
    while ((d & 1) == 0) { s++; d >>= 1; }
    for (i = 0; i < MRMAX; i++){
        if (!millar_rabin(seed_primes[i], s, d, n)) return FALSE;
    }
    return TRUE;
}

#include "bitmap.c"

int main( int argc, char**argv ) {
    U64 start = argc > 1 ? atoll(argv[1]) : 0;
    U64 size  = argc > 2 ? atoll(argv[2]) : 1000*1000*1000;
    U64 i, n;
    char filename[256];
    snprintf(filename, 256, "%llu~%llu.bitmap", start, start+size);
    bitmap * b = bitmap_new(size >> 1 , filename);
    if ((start & 1) == 0) start++;
    for ( i = 0 ; i < size ; i += 2 ) {
	n = start + i;
        if ( is_prime(n) ) {
            bitmap_set( b, i >> 1, 1 );
#ifdef VERBOSE
            printf( "%llu\r", n );
#endif
        }else{
	    bitmap_set( b, i >> 1, 0 );
	}
    }
#ifdef VERBOSE
    printf("\n");
#endif
    return 0;
}