素数10億まで3秒

404 Blog Not Found: C - で素数を数え直したら、範囲10億で10秒切ったお
むむ、以前自分が書いた奴だと、ホットスポットでやってる事はほとんど同じなのに30秒ほどだった。

for (p=2, 3, 5, 7, 11, ...)
   for(i=istart; i<size;i += p*2) pflag[i]=0;

danさんの場合, 1bit でフラグを記憶してるのでメモリが1/8 で済む。そこでメモリアクセスの時間が効いてるんだろう。それならキャッシュに収まる位のブロックに計算を分割しその内側で素数pのループ回せばもっと速くなるかも？と思いやってみた。見事3秒で終わった!
以下コード
danさんのbitmap.cに以下を追加

bitmap *bitmap_block(bitmap *parent, size_t offset, size_t size){
    if (!size) return (bitmap *)NULL;
    if (offset + size > parent->size) return (bitmap *)NULL;

    bitmap *b = (bitmap *)malloc(sizeof(bitmap));
    b->map = &(parent->map[offset>>3]);
    b->size = size;
    b->fd = 0;
    return b;
}

prime.c
単純に計算の順番を変えただけ

// modified by ita 2010/07/29
#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;

// approx cache size
#define BLOCKSIZE ( 1024*1024*4 )

#include "bitmap.c"

#define SIEVEFILE "sieve32.bm"

void sieve_block(bitmap *parent, bitmap *sieve, U64 start0, U64 offset, U64 size_block, U64 pmax)
{
    bitmap *bp = bitmap_block(parent, offset>>1, size_block>>1);
    U64 p, i;
    U64 start = start0 + offset;
    for (p = 3; p < pmax;){
#ifdef VERBOSE
        printf("sieving %llu\r", p);
        fflush(stdout);
#endif
	U64 istart = (( start % p)==0) ? 0 : p - (start % p);
	while( start + istart <= p ) istart += p;
	if ( ((start + istart) &1)==0) istart += p;
        for ( i = istart ; i < size_block ; i += p*2 ) {
            //if (((start + i) & 1) == 0) continue; /* skip even */
            bitmap_set( bp, i>>1, 0 );
        }
        for(p += 2; !bitmap_get(sieve, p>>1); p += 2); /* seek next prime */
    }

    free(bp);
}

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;
    U64 size_orig = size;
    if (size & 0xf)
    {
	size = (size & (~0xfLL))+0x10;
	printf("Size padded: %lld to %lld\n", size_orig, size);
    }
    if (start & 0xf)
    {
	printf("WTF: start & 0xf !=0\n");
	exit(1);
    }
    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;
    U64 block_size = BLOCKSIZE*2;

    for(i = 0; i*block_size < size; i++)
    {
	U64 size1 = block_size;
	if (i*block_size + size1 > size) size1 = size - i*block_size;
        sieve_block(b, sieve, start, i*block_size, size1, pmax);
    }

#ifdef VERBOSE
    // Check number of primes. See http://mathworld.wolfram.com/PrimeCountingFunction.html
    U64 np=0; 
    if (start == 0) np++; // count p=2
    for(i=0; i< size>>1; i++)
    {
	if(bitmap_get(b, i) && (1 + i * 2 <= start + size_orig))
	{
	    np++;
	    //if(np < 200) printf("%lld %lld\n", np, 1+i*2);
	}
    }
    printf("Total %lld primes found\n", np);
#endif

    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;
}