pad.mattdiesel.co.uk

Snippet - Prime Sieve 3 (Friendly)

Prime Sieve 3 (Friendly) (C++)

Prime sieve that's fast and readable now.
Created 2015-07-30 07:33:04.012550 by Matt.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
// I normally compile with
// g++ -std=c++11 -O3 -ffast-math sieve.cc -osieve.exe

#include <iostream>
#include <cmath>

#define DEBUG 1


int applySieve(unsigned int* sieve, const unsigned int n)
{
    unsigned int i, j;
    unsigned int j_max = ((unsigned int)sqrt(n)) >> 1;
    int count = 1; // 2 is a prime.

    // When debug is on, every printevery-th prime is printed.
#if DEBUG
    const int printevery = 10000;
    int x = printevery;
#endif

    for ( i = 1; i < n >> 1; i++ ) {
        if ( !( sieve[i >> 5] & ( 1 << ( i & 31 ) ) ) ) {
            // i*2+1 is a prime

            // Increment prime count
            count++;


            // Decrement print count, when it reaches zero print out the current prime
#if DEBUG
            if ( !--x ) {
                std::cout << ( i << 1 ) + 1 << std::endl;
                x = printevery;
            }
#endif

            // If i is above sqrt(n) then no need to apply the sieve.
            // If count is not required this is where the function can end.
            if (i > j_max) continue;

            // Apply the sieve. Unset any multiples of this prime
            for ( j = ( i << 1 ) + 1 + i ; j < n >> 1; j += ( i << 1 ) + 1 ) {
                sieve[j >> 5] |= 1 << ( j & 31 );
            }
        }
    }

    return count;
}


bool isPrime(const unsigned int* sieve, unsigned int i)
{
    if (i <= 1)
        return false;
    else if (i == 2)
        return true;
    else if (i & 1 == 0) // Sieve uses implicit first bit.
        return false;
    else {
        return !( sieve[i >> 6] & ( 1 << ( (i >> 1) & 31 ) ) );
    }
}


int main( int argc, char const* argv[] )
{
    std::cout << "Prime Sieve 3.3" << std::endl;

    const unsigned int max = 0xFFFFFFFF; // Maximum value for uint32.
    unsigned int* prime = new unsigned int[( max >> 6 ) + 1];

    std::cout << "Memory Allocated." << std::endl;

    int count = applySieve(prime, max);

    std::cout << "Sieve Applied." << std::endl;
    std::cout << "Found " << count << " primes." << std::endl;

    unsigned int test = 0xAB1CDB3; // 10000001st prime
    std::cout << test << " is prime? " << isPrime(prime, test) << std::endl;


    return 0;
}