素数筛法系列之5 完整实现2

// http://acm.nudt.edu.cn/~twcourse/PrimeNumberGeneration.html#a3

#include <math.h>
#include <stdio.h>
#include <memory.h>
#include <time.h>

#define ST             7
#define MOVE           4
#define BLOCKSIZE      (510510 * 4) // 2 * 3 * 5 * 7 * 11 * 13 * 17 = 510510
#define MAXM           ((BLOCKSIZE >> MOVE) + 1)
#define SET_BIT(a, n)  a[n >> 3] |= (1 << (n & 7))

typedef unsigned int uint;
typedef unsigned char uchar;
static uint Prime[6800];
static uchar NumBit0[1 << 8];
static uchar BaseTpl[MAXM * 2];
static const uint MAXN = (-1u / 2);

static int sieve( )
{
    int primes = 1;
    const uint maxp = (1 << 16) + 1000;
    for (uint p = 3; p < maxp; p += 2) {
        if (0 == (BaseTpl[p >> MOVE] & (1 << (p / 2 & 7)))) {
            Prime[primes++] = p;
            for (uint i = p * p / 2; i < maxp / 2; i += p)
                SET_BIT(BaseTpl, i);
        }
    }
    return primes;
}

static void initBit( )
{
    memset(BaseTpl, 0, sizeof(BaseTpl));
    for (int k = 1; k < ST; k++)
    for (int p = Prime[k] / 2; p < sizeof(BaseTpl) * 8; p += Prime[k])
        SET_BIT(BaseTpl, p);

    for (int i = 1; i < sizeof(NumBit0) / sizeof(NumBit0[0]); i++)
        NumBit0[i] = NumBit0[i >> 1] + (i & 1);
    for (int j = 0; j < sizeof(NumBit0) / sizeof(NumBit0[0]); j++)
        NumBit0[j] = 8 - NumBit0[j];
}

static void inline
crossOutFactorMod6(uchar bitarray[], const uint start,
        const uint leng, uint factor)
{
    uint s2, s1 = factor - start % factor;
    s1 += factor - factor * (s1 % 2);
    if (start <= factor)
        s1 = factor * factor;

    const char skip[] = {0, 2, 1, 1, 0, 1};
    const char mrid = ((start + s1) / factor) % 6 - 1;
    s1 = s1 / 2 + skip[mrid] * factor;
    s2 = s1 + skip[mrid + 1] * factor;

    for (factor += 2 * factor; s2 <= leng / 2; ) {
        SET_BIT(bitarray, s1); s1 += factor; //6k + 1
        SET_BIT(bitarray, s2); s2 += factor; //6k + 5
    }
    if (s1 <= leng / 2)
        SET_BIT(bitarray, s1);;
}

static int inline setSieveTpl(uchar bitarray[], uint start, int leng)
{
    const int offset = start % BLOCKSIZE;
    memcpy(bitarray, BaseTpl + (offset >> MOVE), (leng >> MOVE) + 2);
    if (start < 32)
        *((short*)bitarray) = 0x3461; //0x 0011 0100 0110 0001

    bitarray[0] |= (1 << (offset % 16 / 2)) - 1;
    leng += offset % 16;
    *((uint*)bitarray + (leng >> 6)) |= ~((1u << (leng / 2 & 31)) - 1);
    return offset % 16;
}

static int segmentedSieve(uint start, int leng)
{
    uchar bitarray[MAXM + 64];
    const int offset = setSieveTpl(bitarray, start, leng);
    start -= offset, leng += offset;

    const int maxp = (int)sqrt((double)start + leng) + 1;
    for (int i = ST, p = Prime[i]; p < maxp; p = Prime[++i])
        crossOutFactorMod6(bitarray, start, leng, p);

    int primes = 0;
    for (int k = 0; k <= leng >> MOVE; k++)
        primes += NumBit0[bitarray[k]];
    return primes;
}

int countPrimes(uint start, uint stop)
{
    int primes = 0;
    if (start <= 2 && stop >= 2)
        primes = 1;

    const int blocksize = 63 << 14;
    if (start + blocksize > stop)
        return primes + segmentedSieve(start, stop - start + 1);

    primes += segmentedSieve(start, blocksize - start % blocksize);
    for (uint i = start / blocksize + 1; i < stop / blocksize; i++)
        primes += segmentedSieve(i * blocksize, blocksize);

    primes += segmentedSieve(stop - stop % blocksize, stop % blocksize + 1);
    return primes;
}

int main(int argc, char* argv[])
{
    uint start = 0, stop = MAXN;
    clock_t tstart = clock();

    sieve();
    initBit();

    int primeCnt = countPrimes(0, MAXN);
    const char *ouput = "PI(%u, %u) = %d, time use %u ms/n";
    printf(ouput, start, stop, primeCnt, clock() - tstart);

    while (scanf("%u %u", &start, &stop) == 2 && stop >= start) {
        tstart = clock();
        primeCnt = countPrimes(start, stop);
        printf(ouput, start, stop, primeCnt, clock() - tstart);
    }

    return 0;
}
//PI[1, 2147483647] = 105097565, time use 1592 ms
//