segmented sieve wheel base 30

segmented_sieve_wheel_30.cpp

///     This is a implementation of the  wheel segmented sieve with the basis {2,3,5}.
///     This algorithm uses (sqrt(n) *8/30) space.

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

const int64_t L1_CACHE_SIZE = 8192;

const std::vector<int64_t> C1={-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 std::vector<int64_t> C2={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 std::vector<int64_t> RW={-23, -19, -17, -13,  -11,-7 , -1, 1};

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

void Sieve_Wheel_30(int64_t  dim_step, std::vector<char> &Primes)
{
    
    int64_t  m,mmin,i,j;
    int64_t  kmax=(int64_t)std::sqrt(dim_step/bW)+2;
    for (int64_t  k = 1; k  <= kmax; k++)
    {
        for (j = 0; j  < nR; j++)
        {
            if(Primes[k+j*(dim_step+1)])
            {
                for (i = 0; i  < nR; i++)
                {
                    // Cancelling out all the multiples of bW*k+RW[j]   
                    mmin=bW*k*k + k*C1[i*nR+j] + C2[i*nR+j];
                    for ( m =mmin; m <= dim_step; m += bW*k+RW[j])
                        Primes[m+i*(dim_step+1)] = false;
                }
            }
        }
    }
}

void Segmented_Sieve_Wheel_30(int64_t  n)
{
    // counting primes <n     n>bW

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

    if (n>6)
    {
        int64_t  k_end=(n < bW) ? 2 :(int64_t)n/bW+1;
        int64_t  k_sqrt=(int64_t)std::sqrt(k_end/bW)+1;
        int64_t  dim_step,k;   
        int64_t  k_low_2=0;
        dim_step=std::max(L1_CACHE_SIZE-2,k_sqrt);

        std::vector<char> Primes(nR*(dim_step + 1), true);
        Sieve_Wheel_30(dim_step, Primes);

        for (k = 1; k <= std::min(dim_step,k_end-1); k++)
            for (int64_t i = 0; i  < nR; i++)
                if(Primes[k+i*(dim_step+1)])
                    count_p++;
        if (k==k_end && k<=dim_step  && k>0)
            for (int64_t i = 0; i<nR; i++)
                if(Primes[k+i*(dim_step+1)] and RW[i]<(int64_t)(n%bW-bW))
                    count_p++;

        if (k_end>dim_step) 
        {
            int64_t  i,j,kmax,m,mmin,p_j;
            std::vector<char> Primes_t(nR*(dim_step + 1));
 
            for(k_low_2=dim_step;k_low_2<k_end;k_low_2+=dim_step)
            {
                std::fill(Primes_t.begin(), Primes_t.end(), true);

                kmax=std::min((int64_t)std::sqrt((k_low_2+dim_step)/bW)+2,dim_step);
                for(k=1; k<=kmax;k++)
                {
                    for (j = 0; j < nR; j++)
                    {
                        if (Primes[k+j*(dim_step+1)])
                        {
                            for (i = 0; i < nR; i++)
                            {
                                p_j=bW*k+RW[j];
                                mmin=-k_low_2+bW*k*k + k*C1[i*nR+j] + C2[i*nR+j];
                                mmin=(mmin>=0)?mmin:(mmin%p_j+p_j)%p_j;
                                for (m =mmin; m <= dim_step; m += p_j)
                                    Primes_t[m+i*(dim_step+1)] = false;
                            }
                        }
                    }
                }
                for ( k =1+k_low_2; k <=std::min (k_low_2+dim_step,k_end-1); k++)
                    for (i = 0; i < nR; i++)
                        if (Primes_t[k-k_low_2+i*(dim_step+1)])
                            count_p++;
            }
            if (k==k_end && k-k_low_2<=dim_step && k-k_low_2>0)
                for (i = 0; i<nR; i++)
                    if(Primes_t[k-k_low_2+i*(dim_step+1)] && RW[i]<(int64_t)(n%bW-bW))
                        count_p++;
        }
    }

    std::cout << " primes < " << n << ": "<< count_p<< std::endl;
}

int main()
{

    Segmented_Sieve_Wheel_30(100000000);

    return 0;
}