bit segmented sieve wheel base 30

bit_segmented_sieve_wheel_30.cpp

///     This is a implementation of the bit wheel segmented sieve with the basis {2,3,5}.

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <stdint.h>

const int64_t L1_CACHE_SIZE = 32768;

const int64_t bW=30;
const int64_t nR=8;

const int64_t C1[64] = 
{
  -22,  -32,  -28,  -32,  -28,  -8,   -8,   -22,
  -30,  -18,  -30,  -30,  -12,  -30,  -12,  -18,
  -34,  -26,  -16,  -14,  -34,  -26,  -14,  -16,
  -42,  -42,  -18,  -12,  -18,  -18,  -18,  -12,
  -46,  -20,  -34,  -26,  -10,  -14,  -20,  -10,
  -24,  -36,  -36,  -24,  -24,   -6,  -24,   -6,
  -36,  -30,  -24,  -36,  -30,  -24,    0,    0,
  -40,  -38,  -40,  -20,  -22,  -20,   -2,    2
};

const int64_t C2[64] = 
{ 
   0,  9,  7,  9,  7,  1,  1,  0,
   6,  0,  8,  8,  1,  6,  1,  0,
   9,  5,  0,  1,  9,  5,  1,  0,
  15, 15,  1,  0,  3,  3,  1,  0,
  18,  1, 10,  6,  0,  2,  1,  0,
   1, 11, 11,  5,  5,  0,  1,  0,
  10,  7,  4, 10,  7,  4,  0,  0,
  13, 12, 13,  3,  4,  3,  0,  0
};

const int64_t RW[8] = {-23, -19, -17, -13,  -11,-7 , -1, 1};

const int64_t del_bit[8] =
{
  ~(1 << 0),~(1 << 1),~(1 << 2),~(1 << 3),
  ~(1 << 4),~(1 << 5),~(1 << 6),~(1 << 7)
};

const int64_t bit_count[256] =
{
  0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};

void segmented_bit_sieve_wheel_30(uint64_t n)
{

    int64_t k_end = (n < bW) ? (int64_t)2 :(int64_t) (n/(uint64_t)bW+1);
    int64_t k_sqrt = (int64_t) std::sqrt(k_end/bW)+1;
    int64_t segment_size=std::max(L1_CACHE_SIZE-(int64_t)1,k_sqrt); 
    int64_t p_mask_i=(int64_t)4;	

    int64_t  count_p=(n < 3) ? 0 : ((n == 3) ? 1 : ((n <5) ? 2 : 3));

    if (n>6){
        std::vector<uint8_t> Primes(segment_size+(int64_t)1, 0xff);
  
        int64_t  p,m,mmin,i,j,k;
        int64_t kmax = (int64_t) std::sqrt(segment_size/bW)+(int64_t)1;
        for (k = (int64_t)1; k  <= kmax; k++)
        {
            for (j = 0; j  < nR; j++)
            {
                if(Primes[k] & (1 << j))
                {
                    for (i = 0; i<nR; i++)
                    {
                        p=bW*k+RW[j];
                        mmin=bW*k*k + k*C1[i*nR+j] + C2[i*nR+j];
                        for (m =mmin; m <= segment_size && m>=(int64_t)0; m +=p )
                            Primes[m] &= del_bit[i];
                    }
                }
            }
        }
        for (k = (int64_t)1; k <= std::min (segment_size,k_end-(int64_t)1); k++)
            count_p+=bit_count[Primes[k]]; 
        if (k==k_end && k<=segment_size && k>(int64_t)0)
            for (i = 0; i  < nR; i++)
                if(Primes[k]& (1 << i) && RW[i]<(int64_t)(n%bW-bW))
                    count_p++;
 
        if (k_end>segment_size) 
        {
            int64_t  p_j,k_low; 
            std::vector<uint8_t> Segment_t(segment_size+(int64_t)1);
            for (int64_t k_low = segment_size; k_low < k_end; k_low += segment_size)
            {
                std::fill(Segment_t.begin(), Segment_t.end(), 0xff);
                kmax=std::min(segment_size,(int64_t)std::sqrt((k_low+segment_size)/bW)+2);
                for(k=(int64_t)1; k<=kmax;k++)
                {
                    for (j=(int64_t)0; j < nR; j++)
                    {
                        if (Primes[k]& (1 << j))
                        {
                            for (i =(int64_t)0; i < nR; i++)
                            {
                                p_j=bW*k+RW[j];
                                mmin=-k_low+bW*k*k + k*C1[i*nR+j] + C2[i*nR+j];
                                if (mmin<0)
                                    mmin=(mmin%p_j+p_j)%p_j;
                                for (m =mmin; m <= segment_size; m += p_j)
                                    Segment_t[m] &= del_bit[i];
                            }
                        }
                    }
                }
                for ( k =(int64_t)1+k_low; k <=std::min (k_low+segment_size,k_end-(int64_t)1); k++)
                    count_p+=bit_count[Segment_t[k-k_low]];
            }
            if (k==k_end && k-k_low<=segment_size && k-k_low>(int64_t)0)
                for (i = 0; i  < nR; i++)
                    if(Segment_t[k-k_low]& (1 << i) && RW[i]<(int64_t)(n%bW-bW))
                        count_p++;
		}
    }
    std::cout << " primes < " << n << ": "<< count_p<< std::endl;
}

int main()
{
	  segmented_bit_sieve_wheel_30(100000000);
    return 0;
}