ビット配列を使ったエラトステネスの篩を改良

programming c/c++

　ビット配列による素数表よりさらにメモリを節約できるエラトステネスの篩の改良を思いついた。偶数の場合は２以外は素数でないので、素数表にその真偽を保持する必要は無い。なのでインタフェース関数を使うことでメモリ使用量をさらに半分にできる。32bitのunsigned intの配列で表現する場合、Nまでの素数表の作成に必要なメモリはN/32/2になる。
　３以上の奇数の数列を考えると、それぞれを0からインクリメンタルにIDを振るにはnに対してID(n) = n/2-1になる。また、ビット配列のi番目が1かどうか判断するには32bit unsigned intのビット配列tableに対して

idx = i/2-1
item = table[idx / 32] & ((unsigned int)1 << idx % 32)

ビット配列のi番目に0をセットするには

idx = i/2-1
table[idx / 32] &= ~((unsigned int)1 << idx % 32)

となる。
　篩にかけるときは2<=x<=sqrt(N)に対してxの２倍以上の倍数全てを消す。このときに、xの偶数倍はかならず偶数になるが、既に2の倍数として消去されているのでxの奇数倍のときだけ篩にかければよい。xも奇数のときだけ調べればよい。なのでこーなる

for i from 4 to N step 2
  table[i] = 0
for i from 3 to int(sqrt(N)+1) step 2 # 奇数だけ
    for j from i*3 to N step i*2
        table[j] = 0

このアルゴリズムは比較するビット配列版コードで使った。今回の改良されたバージョンではそもそも偶数は素数表に含まれていないので少し異なる。

/*
 * エラトステネスの篩の改良
 * ビット配列に持つのは、3以上の奇数だけ。
 */
#include<cstdio>
#include<climits>
#include<cmath>

#define UINT_BITS (32)
#ifndef N
#define N (1000000)
#endif
#define TBLN (N/UINT_BITS/2)

typedef unsigned int uint;
uint table[TBLN];

bool isprime(int n)
{
    if(n == 2) return true;
    if(n < 2 || (n&1) == 0) return false;
    uint idx = n/2-1;
    uint v = table[idx/UINT_BITS] & ((uint)1 << idx%UINT_BITS);
    return v != 0;
}

void sieve()
{
    for(int i=0; i<=TBLN; i++)
        table[i] = UINT_MAX;
    // エラトステネスの篩
    int u = (int)ceil(sqrt(N)) + 1;
    for(int i=3; i<u; i+=2)
        for(int j=i*3; j<N; j += 2*i){ // j*(偶数)は偶数になるので省く
            int idx = j/2-1;
            table[idx/UINT_BITS] &= ~((uint)1 << idx%UINT_BITS);
        }
}

int main()
{
    sieve();
    return 0;
}

素朴なエラトステネスの篩と比べていろいろやってるので速度もかなり早くなっているはずだと思うので比較してみた。


時間計測はtimeコマンド、環境はdell inspiron mini 9上のvmwareのubuntu 8.04。
コンパイルはg++4.2.4の-O3オプション付きでコンパイルした。
結果：

n 改良ビット配列ver ビット配列ver 素朴ver
1e7 0.728 1.697 6.137
2e7 2.096 3.795 14.971
3e7 3.469 6.053 26.031
4e7 4.774 8.826 31.622
5e7 6.171 11.999 39.622
6e7 7.515 15.148 45.098
7e7 9.172 18.702 58.175
8e7 11.593 20.210 60.906
9e7 12.715 26.519 68.796

グラフ：

メモリが半分ですむだけでなく速くもなった。

以下比較コード：

/*
 * 素朴なエラトステネスの篩
 */
#include<cstdio>
#include<cmath>

#ifndef N
#define N 1000000
#endif

int table[N];

void sieve()
{
    for(int i=0; i<N; i++)
        table[i] = 1;
    table[0] = table[1] = 0;
    int u = (int)ceil(sqrt(N)) + 1;
    for(int i=2; i<u; i++)
        for(int j=i*2; j<N; j+=i)
            table[j] = 0;
}

bool isprime(int n)
{
    return table[n];
}

int main()
{
    sieve();
    return 0;
}

/*
 * ビット配列によるエラトステネスの篩
 */
#include<cstdio>
#include<cmath>
#include<climits>

#define UINT_BITS (32)
#ifndef N
#define N 100
#endif
#define TBLN (N/UINT_BITS)

unsigned int table[N];

void sieve()
{
	for(int i=0; i<=TBLN; i++)
		table[i] = UINT_MAX;
	int u = (int)ceil(sqrt(N)) + 1;
	table[0] &= ~3; // 0,1を0にセット
	for(int i=4; i<N; i+=2)
			table[i/UINT_BITS] &= ~((unsigned int)1 << i%UINT_BITS);
	for(int i=3; i<=u; i+=2)
		for(int j=i*2; j<N; j+=i)
			table[j/UINT_BITS] &= ~((unsigned int)1 << j%UINT_BITS);
}

bool isprime(int n)
{
	unsigned int val = table[n/UINT_BITS] & ((unsigned int)1 << n%UINT_BITS);
	return val != 0;
}

int main()
{
	sieve();
	return 0;
}