Prime numbers using wheel factorial

Question

I have written a program to find prime numbers up to a given limit using wheel factorization of 2, 3 and 5. I would like it to work above \$10^{10}\$, but its still pretty slow.

• What else can I optimize here.
• Any way I can parallelize it?
• I am pretty sure markNonPrimes's algorithm is less efficient than it should be. Is it? I feel that I am brute forcing there: it takes ~250ms when limit is \$10^6\$, ~1500ms when limit is \$10^8\$.

#include <iostream>
#include <vector>
#include <algorithm>
#include <ctime>

typedef uint64_t u64;

const uint8_t ONE = 1;
const u64 F30[] = {1,7,11,13,17,19,23,29};

void markNonPrimes(const u64 start, const u64 num, const u64 lcm, const u64 limit, std::vector<uint8_t> &primes)
{
   for (u64 k = start; (lcm*k) < limit; k++)
   {
      u64 K = lcm*k;
      if ((primes[k] & (ONE << 0)) && (K+F30[0]) % num == 0)
         primes[k] &= ~(ONE << 0);
      if ((primes[k] & (ONE << 1)) && (K+F30[1]) % num == 0)
         primes[k] &= ~(ONE << 1);
      if ((primes[k] & (ONE << 2)) && (K+F30[2]) % num == 0)
         primes[k] &= ~(ONE << 2);
      if ((primes[k] & (ONE << 3)) && (K+F30[3]) % num == 0)
         primes[k] &= ~(ONE << 3);
      if ((primes[k] & (ONE << 4)) && (K+F30[4]) % num == 0)
         primes[k] &= ~(ONE << 4);
      if ((primes[k] & (ONE << 5)) && (K+F30[5]) % num == 0)
         primes[k] &= ~(ONE << 5);
      if ((primes[k] & (ONE << 6)) && (K+F30[6]) % num == 0)
         primes[k] &= ~(ONE << 6);
      if ((primes[k] & (ONE << 7)) && (K+F30[7]) % num == 0)
         primes[k] &= ~(ONE << 7);
   }
}

void getPrimesUpto(u64 limit)
{
   const u64 lcm = 30;
   std::vector<uint8_t> primes(limit*0.26 + 1,0xff);

   u64 k = 1;
   u64 n = 0;
   u64 count=3-1;

   markNonPrimes(k, F30[1], lcm, limit, primes);
   markNonPrimes(k, F30[2], lcm, limit, primes);
   markNonPrimes(k, F30[3], lcm, limit, primes);
   markNonPrimes(k, F30[4], lcm, limit, primes);
   markNonPrimes(k, F30[5], lcm, limit, primes);
   markNonPrimes(k, F30[6], lcm, limit, primes);
   markNonPrimes(k, F30[7], lcm, limit, primes);

   while(lcm*k*lcm*k<limit)
   {
      const u64 K = lcm*k;
      if(primes[k] & (ONE<<0))
         markNonPrimes(k+1, K+F30[0], lcm, limit, primes);

      if(primes[k] & (ONE<<1))
         markNonPrimes(k+1, K+F30[1], lcm, limit, primes);

      if(primes[k] & (ONE<<2))
         markNonPrimes(k+1, K+F30[2], lcm, limit, primes);

      if(primes[k] & (ONE<<3))
         markNonPrimes(k+1, K+F30[3], lcm, limit, primes);

      if(primes[k] & (ONE<<4))
         markNonPrimes(k+1, K+F30[4], lcm, limit, primes);

      if(primes[k] & (ONE<<5))
         markNonPrimes(k+1, K+F30[5], lcm, limit, primes);

      if(primes[k] & (ONE<<6))
         markNonPrimes(k+1, K+F30[6], lcm, limit, primes);

      if(primes[k] & (ONE<<7))
         markNonPrimes(k+1, K+F30[7], lcm, limit, primes);  

      k++;
   }

   for(u64 k=0; lcm*k < limit; k++ )
   {
      for(u64 i =0; i<8; i++)
         if(primes[k] & ONE<< i && lcm*k+F30[i] < limit) {count++;/*std::cout << lcm*k+F30[i] << " ";*/}
   }

   std::cout  << "\n======\n"<< count << "\n";
}

int main()
{
   clock_t start = clock();
   getPrimesUpto(100);
   std::cout << "\nTime: "<< clock() - start << std::endl;
}

Answer 1

1. for loops

   You are doing the same thing over and over again in many places in the code.

   For example,

   for (u64 k = start; (lcm*k) < limit; k++){
       u64 K = lcm*k;
       if ((primes[k] & (ONE << 0)) && (K+F30[0]) % num == 0)
           primes[k] &= ~(ONE << 0);
           ......
           if ((primes[k] & (ONE << 7)) && (K+F30[7]) % num == 0)
               primes[k] &= ~(ONE << 7);
   }
   

   could become

   for(u64 k = start; lcm*k < limit; k++)
   {
       u64 K = lcm*k;
       for(int i = 0; i < 8; i++)
       {
           if ((primes[k] & (ONE << i)) && (K+F30[i]) % num == 0)
               primes[k] &= ~(ONE << i);
       }
   }
   

   There are many other places you could do this as well. Also, the while loop near the bottom could be transformed into a for loop, as you are basically doing a for loop with while.

2. u64 vs int

   Don't use a u64 when you don't need them; ints save on storage. You can replace the bottom two loops with regular ints.

3. Comments

   Your code is somewhat confusing when you look at it on first glance. Maybe you should add some comments for the reader to see what is happening in your MarkNonPrimes function. Also, in the variable count (why declare it to 3-1 instead of 2?), and how you got that magic constant 0.26 in your vector declaration. The naming of F30 is pretty cryptic too.

4. Conventions

   If you are going to use a convention then stick with it. Don't use ONE << 0 then ONE<<0 later.

5. Global variables

   See here. Pass the variables as a parameter instead of declaring them globally, because as the project grows in size, then it will become harder to not have an ambiguous match for your global names. However, the constant ONE is fine, because not many variables will have the name one, but you should rename it u8_ONE instead because you may be declaring a u16_ONE or something of the sort.

Here is the code with all the modifications (except comments and naming):

#include <iostream>
#include <vector>
#include <algorithm>
#include <ctime>

typedef uint64_t u64;
const uint8_t u8_ONE = 1;    
void markNonPrimes(const u64 start, const u64 num, const u64 lcm, const u64 limit, std::vector<uint8_t> &primes, u64 F30[])
{
    for(u64 k = start; lcm*k < limit; k++)
    {
        u64 K = lcm*k;
        for(int i = 0; i < 8; i++)
        {
            if ((primes[k] & (u8_ONE << i)) && (K + F30[i]) % num == 0)
                primes[k] &= ~(u8_ONE << i);
        }
    }
}

void getPrimesUpto(u64 limit)
{
   const u64 F30[] = {1,7,11,13,17,19,23,29};

   const u64 lcm = 30;
   std::vector<uint8_t> primes(limit*0.26 + 1,0xff);
   u64 n = 0;
   u64 count = 2;
   for(int i = 1; i < 8; i++)
   {
       markNonPrimes(1, F30[i], lcm, limit, primes, F30);
   }
   for(int k = 1; lcm * k * lcm * k < limit; k++)
   {
      const u64 K = lcm * k;
      for(int i = 0; i < 8; i++)
          if(primes[k] & (u8_ONE << i))
              markNonPrimes(k+1, K+F30[i], lcm, limit, primes, F30);
   }

   for(int k=0; lcm*k < limit; k++)
   {
      for(int i = 0; i < 8; i++)
         if(primes[k] & u8_ONE << i && lcm * k + F30[i] < limit)  
         {
             count++;
             /*std::cout << lcm*k+F30[i] << " ";*/
         }
   }

   std::cout  << "\n======\n"<< count << "\n";
}

int main()
{
   clock_t start = clock();
   getPrimesUpto(100);
   std::cout << "\nTime: "<< clock() - start << std::endl;
}