Reduce Prime Sieve Memory Consumption 2

Question

I've made a prime number generator (for Project Euler). It uses Euler's Sieve (a modified Sieve of Eratosthenes), with a mod 30 step. I'd like to reduce the memory consumption to 4/15 what it currently is by keeping a boolean array only for the possibly prime remainders of 30. I can't get it to work and also fear that this will slow down the program.

I've used

n[f/30*8+[zeroes, except 1,2,3,4,5,6,7 at 7,11,13,17,19,23,29][f%30]]

which seemed to not filter out any composite numbers. How can I make this work and what other optimizations (aside from increasing the mod) can you suggest? I need the primes up to about two billion.

//pes30.cxx
#include <vector>
#include <algorithm>
#ifndef _pes30_cxx_
#define _pes30_cxx_

typedef unsigned long long big;

const int offsets[]={6,4,2,4,2,4,6,2};

//fill a given vector with all primes under some value:
void sievePrimes(big max, std::vector<big> &p){
    big multiple;
    bool n[max];//array of whether or not each number is prime (30/8 too big!)
    p={2,3,5};//because the sieve skips all multiples of 2,3, and 5, start with them.
    for(big i=0; i<max; ++i)//initialize the array
        n[i]=true;
    //for every number marked prime less than max, mark its multiples with
    //    every number still marked prime over it as composite.
    for(big i=7, step=1; i<=max; i+=offsets[step], ++step==8?step=0:0){
        if(!n[i])//if i is not prime
            continue;
        p.push_back(i);//add i to the list of primes
        //finds every multiple of i and a (still marked) prime greater than i
        for(big j=i, step2=step; j<=max/i; j+=offsets[step2], ++step2==8?step2=0:0){
            if(!n[j])//skip nonprimes
                continue;
            multiple=j*i;//begin at i^2
            do{
                n[multiple]=false;
            }while((multiple*=j) <= max);
        }
    }
}
//test if a number is prime by searching a given list of primes for it:
inline bool isPrime(big n, std::vector<big> p){
    return std::binary_search(p.begin(), p.end(), n);
}
#endif

Note that this is a re-post with ... proper formatting, and while I appreciate critiques of my coding style, I would also appreciate some advice about improving the array as well (storing only numbers n with n mod 30 prime or 1, but not 2, 3, or 5).

Answer 1

You've posted some code, which shows what you told the computer to do. And it even has some comments. Honestly, that's a lot better than the average question.

Then you realize that we aren't going to understand this blob of code, so you use phrases like "I made" and "It uses Euler's Sieve" and "I need the primes up to about two billion". Look, that's super-useful information, and I'm glad you told me that, but I'm going to suggest that in the future you put information about who wrote the program, the name of the algorithm you use, etc. into comments embedded in the program:

// pes30.cxx
// prime number generator for Project Euler
// Print all the primes up to 2x10^9.
// 2014-05-08: started by hacatu
// implementation of Euler's sieve (a modified Sieve of Eratosthenes)
// http://en.wikipedia.org/wiki/Euler%27s_sieve
// using wheel factorization with a wheel of size 30 to reduce memory requirements
// http://en.wikipedia.org/wiki/wheel_factorization

So is there some reason not to use http://primesieve.org/ ? Have you gotten it working without the wheel factorization optimization? Have you tried looking at only odd prime candidates (wheel of size 2) or a small wheel of size 6 or 12 ?

I am mystified by the 3 lines

bool n[max];//array of whether or not each number is prime (30/8 too big!)
for(big i=0; i<max; ++i)//initialize the array
    n[i]=true;

The "max" is the largest number you want to check for possibly being prime during this run of the program, right? So that line eats a bunch of memory -- one bool per each number from 0 to max, perhaps 2 million bool values.

I thought you were going to use wheel factorization to reduce memory requirements?

I think you want to replace those 3 lines with something more like the following 2 lines of code:

big blocks_of_30_values = (max/30)+1;
// For example, when max is 38, max/30 gives 1, but we need 2 blocks: 0..29 and 30..59.
// In each block of 30 values under consideration
// (i.e, the block 30..59, the block 60..89, the block 90..119, etc.),
// after using the wheel factorization to eliminate multiples of 2 and 3 and 5,
// there are only 8 potential candidates left in each block.
// So store each block in a single unsigned char:
// and initialize the array to all 1 bits (1=candidate prime, 0=definitely composite)
std::vector<unsigned char> candidatePrime( blocks_of_30_values, 0xFF );

So when max is 2 billion, we only allocate (max/30) = under 70 million chars, or roughly 533 million bits.

It might make your program easier to read if you have named functions for converting a bit in the array to the integer value it represents and vice versa:

// the bit at the n'th bit of the m'th byte of the candidatePrime array
// represents what value?
big the_value( big block_number, int bit_number ){
    const int offset[] = { 1, 7, 11, 13, 17, 19, 23, 29 };
    int this_offset = offset[bit_number];
    return( (30*block_number) + this_offset );
}
// the given integer value n corresponds to
// what (block_number, bit_number) of the candidatePrime array?
// assumes that n is *not* a multiple of 2, 3, or 5.
big block_number( big integer_value ){
    return (integer_value/30);
}
int bit_number( big integer_value ){
    // as long as the integer_value is not a multiple of 2, 3, or 5,
    // this function should never return '9'.
    // translate a value 0..29 into the corresponding bit number 0..7
    const int bit_offset[] = {
        9, 0, 9,  9, 9, 9,
        9, 1, 9,  9, 9, 2,
        9, 3, 9,  9, 9, 4,
        9, 5, 9,  9, 9, 6
        9, 7, 9,  9, 9, 7
    };
    return bit_offset[ integer_value % 30 ];
}

Answer 2

I have fixed the vector to use one thirtieth the memory by using offsets and also by using a vector<boolean>, which is optimized to use every bit. This caused a 50% slowdown, but it is worth it for high max values (because the max value can now be 30 times as great), and I'm sure I can reduce this.

I did not use the primesieve.org because I did not know about it. It is far better than my program. It generates the primes under 1 billion in 0.169s, my program takes 18.539s.

Here is my updated program, now with a sensibly sized vector (and even more overworked for loop increments):

// pes30.cxx
// prime number generator for Project Euler
// Print all the primes up to 2x10^9.
// 2014-05-08: started by hacatu
// implementation of Euler's sieve (a modified Sieve of Eratosthenes)
// http://en.wikipedia.org/wiki/Euler%27s_sieve
// using wheel factorization with a wheel of size 30 to reduce memory requirements
// http://en.wikipedia.org/wiki/wheel_factorization
#include <vector>
#include <algorithm>
#include <cstdint>
#include <cmath>
#ifndef _pes30_cxx_
#define _pes30_cxx_
using namespace std;
const int offsets[]={6,4,2,4,2,4,6,2};//steps in the wheel 1,7,11,13,17,19,23,29,1...
const int bool_offsets[]={//used to find a number's spot on the wheel, for accessing its primality in the vector
9, 0, 9, 9, 9, 9,//the nines should never be accessed
9, 1, 9, 9, 9, 2,
9, 3, 9, 9, 9, 4,
9, 5, 9, 9, 9, 6,
9, 9, 9, 9, 9, 7
};
//create a vector with all primes under some value:
inline vector<uint64_t> sievePrimes(uint64_t max){
    uint64_t multiple;
    vector<bool> n(max/30*8+8, true);//vector of whether every number relatively prime to 30 is still a candidate prime.
    //8 bits are needed for every 30 numbers (since 8 are relatively prime with 30), so max/30*8, plus 8 because max/30 rounds down.
    vector<uint64_t> primes={2,3,5};//because the sieve skips all multiples of 2,3, and 5, start with them.
    //for every number marked prime less than max, mark its multiples with
    //every number still marked prime over it as composite.
    int r=sqrt(max);
    for(uint64_t i=1, p=7, step=1; p <= max; ++i, p += offsets[step], ++step == 8 ? step=0:0){
        if(!n[i])//if p is not prime (using i for the index holding the primality of p, to avoid computing the index for every number)
            continue;
        primes.push_back(p);//add p to the list of primes
        //finds every multiple of i and a (still marked) prime greater than i
        for(uint64_t j=i, p2=p, step2=step; p2 <= max/p; ++j, p2 += offsets[step2], ++step2 == 8 ? step2=0:0){
            if(!n[j])//skip nonprimes
                continue;
            multiple=p*p2;
            do{
                n[multiple/30*8 + bool_offsets[ multiple%30 ]]=false;
            }while((multiple*=p) <= max);
        }
    }
    return primes;
}
//test if a number is prime by searching a given list of primes for it:
inline bool isPrime(uint64_t n, vector<uint64_t> primes){
    return std::binary_search(primes.begin(), primes.end(), n);
}
#endif