35msくらいで100万番目の素数までのエラトステネスの篩

main.d

/**
# 35msくらいで100万番目の素数を求める
$ ldmd2 -O -release -inline foo.d
15485863
34
参考：http://d.hatena.ne.jp/k3_kaimu/20120806/1344253324
*/
import std.stdio;
import std.math;
import std.datetime : StopWatch;
import core.bitop   : bts, popcnt, bt;

enum long M = 10_000_000;
enum long N = M * 18;
enum Leng = N/3 + (N%6 == 2);
enum NB = size_t.sizeof * 8;

void main(){
    StopWatch sw;
    sw.start;
    
    int k, k2;
    size_t[] buf = new size_t[Leng / NB];
    size_t* p = buf.ptr;
    
    int cnt = 1, idx;
    
    foreach(i; 1..cast(uint)((N^^0.5)/3 + 1)){
        k = (3 * i + 1) | 1;
        k2 = k << 1;
        
        for(int j = k*k/3; j < Leng; j += k2)
            bts(p, j);
        for(int j = k*(k-((i&1)<<1)+4)/3; j < Leng; j += k2)
            bts(p, j);
    }
    
    
    long result;
    LOOP:
    foreach(e; p[0..Leng/NB]){
        
        cnt += NB - popcnt(e);
        
        if(cnt >= M){
            cnt -= NB - popcnt(e);
            idx *= NB;
            
            for(;;++idx){
                cnt += !bt(p, idx);
                
                if(cnt == M){
                    result = (3 * idx + 1)|1;
                    break LOOP;
                }
            }
        }
        ++idx;
    }
    sw.stop();

    writeln(result);
    writeln(sw.peek.msecs);
}