- #include <math.h>
- #include <stdio.h>
- #include <memory.h>
- #include <time.h>
- #define ST 7
- #define MOVE 4
- #define BLOCKSIZE (510510 * 8) // 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 short ushort;
- typedef unsigned char uchar;
- static ushort NumBit0[1 << 16];
- static uint Prime[6800];
- 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] = 16 - 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) = 0x3491; //0x 0011 0100 1001 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 countBit0(ushort* bitarray, int wsize)
- {
- int primes = 0;
- while (wsize-- >= 0)
- primes += NumBit0[*bitarray++];
- return primes;
- }
- 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);
- return countBit0((ushort*)bitarray, leng >> MOVE >> 1);
- }
- int countPrimes(uint start, uint stop)
- {
- int primes = 0;
- if (start <= 2 && stop >= 2)
- primes = 1;
- const int segsize = 63 << 14;
- if (start + segsize > stop)
- return primes + segmentedSieve(start, stop - start + 1);
- primes += segmentedSieve(start, segsize - start % segsize);
- for (uint i = start / segsize + 1; i < stop / segsize; i++)
- primes += segmentedSieve(i * segsize, segsize);
- primes += segmentedSieve(stop - stop % segsize, stop % segsize + 1);
- return primes;
- }
- int main(int argc, char* argv[])
- {
- uint start = 0, stop = MAXN;
- clock_t tstart = clock();
- sieve();
- initBit();
- int primeCnt = countPrimes(0, stop);
- 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
- //
- //该片段来自于http://www.codesnippet.cn/detail/060820134992.html
来源: http://www.codesnippet.cn/detail/060820134992.html