Segmented Sieve Wheel base 6

Segmented_Sieve_Wheel_6.cpp

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

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

const int64_t L1_CACHE_SIZE = 32768;

void Sieve_Wheel_6(int64_t  dim_step, std::vector<char> &Primes5mod6, std::vector<char> &Primes1mod6)
{
    // Get all the primes >3 and <= 6*dim_step+1 in vectors Primes5mod6 and Primes1mod6
    // 6*i-1 or 6*i+1 are prime if [i] element is true in Primes5mod6 or Primes1mod6
    // the sieve only considers numbers +-1 mod 6 a kind of Wheel factorization with the basis 2,3
    
    int64_t  m,mmin;
    int64_t  imax=(int64_t)std::sqrt(dim_step/6)+2;
    for (int64_t  i = 1; i  <= imax; i++)
    {
        if (Primes5mod6[i])
        {
            // Cancelling out all the multiples of 6*i-1
            mmin=6*i*i;
            for ( m =mmin; m <= dim_step; m += 6*i-1)
            {
                Primes5mod6[m] = false;
                Primes1mod6[m-2*i] = false;
            }
            for (; m-2*i <= dim_step; m += 6*i-1)
                Primes1mod6[m-2*i] = false;
        }
        if (Primes1mod6[i])
        {
            // Cancelling out all the multiples of 6*i+1
            mmin=6*i*i;
            for ( m =mmin; m+2*i <= dim_step; m += 6*i+1)
            {
                Primes5mod6[m] = false;
                Primes1mod6[m+2*i] = false;
            }
            for (; m <= dim_step; m += 6*i+1)
                Primes5mod6[m] = false;
       }
    }
}

void Segmented_Sieve_Wheel_6(int64_t  n)
{
    // counting primes <n

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

    if (n>=6)
    {
        int64_t  k_end=(int64_t)n/6;
        int64_t  k_sqrt=(int64_t)std::sqrt(k_end/6)+1;
        int64_t  dim_step=(int64_t)exp2(10)-1;   
  
        //For small numbers  dim_step reduced to 2^10-1  
        if (k_sqrt>dim_step) 
            dim_step=std::max(L1_CACHE_SIZE-2,k_sqrt);

        // Get all the primes >3 <= 6*dim_step+1 
        std::vector<char> Primes5mod6(dim_step + 1, true);
        std::vector<char> Primes1mod6(dim_step + 1, true);
        Sieve_Wheel_6(dim_step, Primes5mod6, Primes1mod6);

        // count primes in range [5, 6*dim_step+1] 
        for (int64_t k = 1; k <= std::min(dim_step,k_end); k++)
        {
            if (Primes5mod6[k])
                count_p++;
            if (Primes1mod6[k])
                count_p++;
        }

        if (k_end>dim_step) 
        {
            int64_t  i,imax,k,k_low_2,m,m1,mmin0,mmin,mmin1,p_i;
            std::vector<char> Primes5mod6_t(dim_step + 1);
            std::vector<char> Primes1mod6_t(dim_step + 1);  
            for(k_low_2=dim_step;k_low_2<k_end;k_low_2+=dim_step)
            {
                std::fill(Primes5mod6_t.begin(), Primes5mod6_t.end(), true);
                std::fill(Primes1mod6_t.begin(), Primes1mod6_t.end(), true);

                imax=std::min((int64_t)std::sqrt((k_low_2+dim_step)/6)+2,dim_step);
                for(i=1; i<=imax;i++)
                {
                    if (Primes5mod6[i])
                    {
                        // Cancelling out all the multiples of p_i=6*i-1 in Primes5mod6_t and Primes1mod6_t
                        p_i=6*i-1;
                        mmin0=-k_low_2+6*i*i;
                        mmin1=(mmin0-2*i>=0)?mmin0-2*i:((-k_low_2-i)%p_i+p_i)%p_i;
                        mmin=(mmin0>=0)?mmin0:(mmin1+2*i)%p_i;
                        for (m =mmin,m1 =mmin1; (m <= dim_step)&&(m1 <= dim_step); m += p_i,m1 += p_i)
                        {
                            Primes5mod6_t[m] = false;
                            Primes1mod6_t[m1] = false;
                        }
                        for (; m <= dim_step; m += p_i)
                            Primes5mod6_t[m] = false;
                        for (; m1 <= dim_step; m1 += p_i)
                            Primes1mod6_t[m1] = false;
                    }
                    if (Primes1mod6[i])
                    {
                        // Cancelling out all the multiples of p_i=6*i+1 in Primes5mod6_t and Primes1mod6_t
                        p_i=6*i+1;
                        mmin0=-k_low_2+6*i*i;
                        mmin=(mmin0>=0)?mmin0:((-k_low_2-i)%p_i+p_i)%p_i;
                        mmin1=(mmin0+2*i>=0)?mmin0+2*i:(mmin+2*i)%p_i;
                        for (m =mmin,m1 =mmin1; (m <= dim_step)&&(m1 <= dim_step); m += p_i,m1 += p_i)
                        {
                            Primes5mod6_t[m] = false;
                            Primes1mod6_t[m1] = false;
                        }
                        for (; m <= dim_step; m += p_i)
                            Primes5mod6_t[m] = false;
                        for (; m1 <= dim_step; m1 += p_i)
                            Primes1mod6_t[m1] = false;
                    }
                }
                for ( k =1+k_low_2; k <=std::min (k_low_2+dim_step,k_end); k++)
                {
                    if (Primes5mod6_t[k-k_low_2])
                        count_p++;
                    if (Primes1mod6_t[k-k_low_2])
                        count_p++;
                }
            }
            if (k==k_end+1 and k_end-k_low_2<=dim_step)
                if(n%6<=1 and Primes1mod6_t[k_end-k_low_2])
                	count_p--;
        }
        else
        {
             if(n%6<=1 and Primes1mod6[k_end])
                count_p--;
        }
    }
    std::cout << " primes < " << n << ": "<< count_p<< std::endl;
}

int main()
{

    Segmented_Sieve_Wheel_6(100000000);

    return 0;
}