5
\$\begingroup\$

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).

\$\endgroup\$
0

2 Answers 2

4
\$\begingroup\$

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 ];
}
\$\endgroup\$
2
\$\begingroup\$

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
\$\endgroup\$
7
  • \$\begingroup\$ vector<bool> was 50% slower? Did you compile with optimization enabled? Or probably it uses 64-bit chunks, so every RMW to update a bit is an RMW of the whole 64-bit integer. So crossing off small primes creates long chains of store/reload latency as many indices map to the same 64-bit chunk. (8-bit chunks have an advantage for this, but big downsides for other cases like find_if. See Sieve of Eratosthenes in x86 assembly - I actually got a speedup using 16-bit operand-size in hand-written code in ways a compiler might not. \$\endgroup\$ Aug 25, 2022 at 1:19
  • \$\begingroup\$ That was on a non-segmented non-wheel sieve, trying to keep things pretty simple and mostly micro-optimize; with more work to figure out where to mark off each bit, it would be a bit less sensitive to those details. \$\endgroup\$ Aug 25, 2022 at 1:23
  • 1
    \$\begingroup\$ I don't remember if I compiled with optimizations or not, but yes I expect that the "bitvector" bundling 64 bits together in a read-heavy scenario has to do with it. This is OLD code, and it has MAJOR problems: (1) does not correctly check for overflow (2) the wheel should be implemented OUTSIDE the main loop, to remove branching from the inner loop. This is the key insight that really fast implementations like primesieve use (3) this isn't actually segmented (4) std::vector doesn't inform the compiler that these containers never grow (5) binary_search in isPrime is bizarre \$\endgroup\$
    – hacatu
    Aug 25, 2022 at 7:07
  • 1
    \$\begingroup\$ xD, I hadn't read the code in detail. Yeah, making a vector<uint64_t> primes is self-defeating vs. just keeping the bitmap to answer isPrime queries from. Even for answering next_prime queries to find a prime above a certain start point, linear search of a bitmap is good; primes are dense enough. Especially if the library has a good specialization of std::find for vector<bool> bitmaps, checking whole chunks at a time (isocpp.org/blog/2012/11/on-vectorbool). \$\endgroup\$ Aug 25, 2022 at 11:40
  • 1
    \$\begingroup\$ The isPrime function was secondary; I still mostly wanted a list of primes, but using binary search on that list to check for primes is terrible. For one thing, interpolation search could be used. For another, even on the whole 64 bit range, storing the bitvector is better for isPrime since there's only n/log(n) primes up to n plus a bitvector is "O(1)" addressable. Using the MIller Rabin test with a few fixed bases is probably the best way to do prime testing in the 64 bit range though, which is mostly what I was thinking of when I said "binary_search in isPrime is bizarre". \$\endgroup\$
    – hacatu
    Aug 25, 2022 at 19:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.