Segmented wheel sieve

Double_Segmented_Wheel_Sieve.cpp

///     This is a implementation of the  wheel segmented sieve that 
///     uses two vectors in order to have the wheel with the basis {2,3}.
///     This sieve 1/3 of the memory is used compared to the traditional sieve.
///     For large numbers a double segmentation is used.

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

const int64_t L1_CACHE_SIZE = 32768;

void SieveWheel6(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;
    int64_t  mmin=6; //mmin=6*i*i with i=1
    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
            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
            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;
       }
       mmin+=12*i+6; //mmin=6*(i+1)*(i+1)=6*i*i+12*i+6=mmin+12*i+6
    }
}

void StepSegmentedSieve(int64_t  k_low,int64_t  dim_step,int64_t  k_low_2,int64_t  dim_step_2, std::vector<char> &Primes5mod6, std::vector<char> &Primes1mod6, std::vector<char> &Primes5mod6_t, std::vector<char> &Primes1mod6_t)
{
    int64_t  imax=(int64_t )std::sqrt((k_low_2+dim_step_2)/6)+2;
    int64_t  m,m1,mmin,mmin1,p_i;
    int64_t  mmin0=6-k_low_2+12*k_low+6*k_low*k_low;
    imax=std::min(imax,dim_step);
    for(int64_t i=1; i<=imax;i++)
    {
        if (Primes5mod6[i])
        {
            // Cancelling out all the multiples of p_i=6*(i+k_low)-1
            p_i=6*(i+k_low)-1;
	    mmin1=(mmin0-2*(i+k_low)>=0)?mmin0-2*(i+k_low):((-k_low_2-i-k_low)%p_i+p_i)%p_i;
	    mmin=((mmin0)>=0)?mmin0:(mmin1+2*(i+k_low))%p_i;
            for (m =mmin,m1 =mmin1; (m <= dim_step_2)&&(m1 <= dim_step_2); m += p_i,m1 += p_i)
            {
                Primes5mod6_t[m] = false;
	        Primes1mod6_t[m1] = false;
            }
	    for (; m <= dim_step_2; m += p_i)
                Primes5mod6_t[m] = false;
	    for (; m1 <= dim_step_2; m1 += p_i)
                Primes1mod6_t[m1] = false;
        }
        if (Primes1mod6[i])
        {
            // Cancelling out all the multiples of p_i=6*(i+k_low)+1
            p_i=6*(i+k_low)+1;
	    mmin=(mmin0>=0)?mmin0:((-k_low_2-i-k_low)%p_i+p_i)%p_i;
	    mmin1=((mmin0+2*(i+k_low))>=0)?mmin0+2*(i+k_low):(mmin+2*(i+k_low))%p_i;
            for (m =mmin,m1 =mmin1; (m <= dim_step_2)&&(m1 <= dim_step_2); m += p_i,m1 += p_i)
            {
                Primes5mod6_t[m] = false;
	        Primes1mod6_t[m1] = false;
            }
	    for (; m <= dim_step_2; m += p_i)
                Primes5mod6_t[m] = false;
	    for (; m1 <= dim_step_2; m1 += p_i)
                Primes1mod6_t[m1] = false;
        }
        mmin0+=12*i+6+12*k_low;
    }
}

void DoubleSegmentedSieve( int64_t  k_start, int64_t  k_end,int64_t dim_step_1,int64_t dim_step_2, std::vector<uint64_t> &ans)
{
    
    int64_t  sqrt_max=(int64_t)std::sqrt(k_end/6)+2;
    int64_t  dim_step=(int64_t)exp2(10)-1;   
    int64_t  k_low,k_low_2;
    int64_t  num_step=0;

    //For small numbers  dim_step_1 reduced to 2^10  
    if (sqrt_max>dim_step) 
        dim_step=dim_step_1;
        
    // Get all the primes >3 <= 6*dim_step+1 for double segmented is always 6*dim_step+1 >sqrt(sqrt(6*k_end+1))
    std::vector<char> Primes5mod6(dim_step + 1, true);
    std::vector<char> Primes1mod6(dim_step + 1, true);
    SieveWheel6(dim_step, Primes5mod6, Primes1mod6);

    // Append  in ans primes in range [6*k_start, 6*dim_step+1] if k_start<dim_step
    int64_t  k_min=k_start+1;
    if (k_start==0)
    {
        ans.push_back(2);
        ans.push_back(3);
        ans.push_back(5);
        ans.push_back(7);
        k_min++;
    }
    if (k_start<=dim_step&&(k_min==k_start+1))
        if (Primes1mod6[k_start])
            ans.push_back((uint64_t)6*k_start+1); 
    for (int64_t k = k_min; k <= std::min(dim_step,k_end); k++)
    {
        if (Primes5mod6[k])
            ans.push_back((uint64_t)6*k-1);
        if (Primes1mod6[k])
            ans.push_back((uint64_t)6*k+1);
    }

    if (sqrt_max>dim_step_1) //Double Segmented
    {
        int  ck_ans=0;
        num_step=(int64_t)sqrt_max/dim_step;
        if (sqrt_max>num_step*dim_step)
            num_step+=1;
        
        //prime temp vectors in a single step [low,low+dim_step_2] in the interval [sqrt_max+c+1,k_end]
        std::vector<char> Primes5mod6_f(dim_step_2 + 1, true);
        std::vector<char> Primes1mod6_f(dim_step_2 + 1, true);
        
        //prime temp vectors in a single step [j*dim_step,(j+1)*dim_step] in the interval [0,sqrt_max+c] with (sqrt_max+c)%dim_step=0
        std::vector<char> Primes5mod6_t(dim_step + 1);
        std::vector<char> Primes1mod6_t(dim_step + 1); 

        for(k_low_2=std::max(num_step*dim_step,k_start-1);k_low_2<k_end;k_low_2+=dim_step_2)
        {
            //removes the multiples of the primes of Primes5mod6 Primes1mod6           
            StepSegmentedSieve(0,dim_step,k_low_2, dim_step_2, Primes5mod6, Primes1mod6,Primes5mod6_f, Primes1mod6_f);

            for(k_low=dim_step;k_low<sqrt_max;k_low+=dim_step)
            {
                std::fill(Primes5mod6_t.begin(), Primes5mod6_t.end(), true);
	        std::fill(Primes1mod6_t.begin(), Primes1mod6_t.end(), true);
	        
                //find the primes of Primes5mod6_t Primes1mod6_t
                StepSegmentedSieve(0,dim_step,k_low, dim_step,Primes5mod6, Primes1mod6,Primes5mod6_t, Primes1mod6_t);
                
                // if k_start<sqrt_max+c get in ans primes in the interval [k_start,sqrt_max+c]
                if (ck_ans==0)
                {
                    for (int64_t k =std::max(1+k_low,k_start-1); k <=std::min (k_low+dim_step,k_end); k++)
                    {
                        if (Primes5mod6_t[k-k_low])
                            ans.push_back((uint64_t)6*k-1);
                        if (Primes1mod6_t[k-k_low])
                           ans.push_back((uint64_t)6*k+1);
                    }
                }

                //removes the multiples of the primes of Primes5mod6_t Primes1mod6_t
                StepSegmentedSieve(k_low,dim_step,k_low_2, dim_step_2, Primes5mod6_t, Primes1mod6_t,Primes5mod6_f, Primes1mod6_f);
            }
            ck_ans=1;
            //get in ans  the primes of Primes5mod6_f Primes1mod6_f
            for (int64_t k =1+k_low_2; k <=std::min (k_low_2+dim_step_2,k_end); k++)
            {
                 if (Primes5mod6_f[k-k_low_2])
                     ans.push_back((uint64_t)6*k-1);
                 if (Primes1mod6_f[k-k_low_2])
                     ans.push_back((uint64_t)6*k+1);
            }
            std::fill(Primes5mod6_f.begin(), Primes5mod6_f.end(), true);
	    std::fill(Primes1mod6_f.begin(), Primes1mod6_f.end(), true);
        }
    }
    else if (k_end>dim_step) //Single Segmented    or without segmentation  for small numbers 
    {
        std::vector<char> Primes5mod6_t(dim_step + 1);
        std::vector<char> Primes1mod6_t(dim_step + 1);  
        for(k_low_2=std::max(dim_step,k_start-1);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);
	    StepSegmentedSieve(0,dim_step,k_low_2, dim_step, Primes5mod6, Primes1mod6,Primes5mod6_t, Primes1mod6_t);
	    for (int64_t k =1+k_low_2; k <=std::min (k_low_2+dim_step,k_end); k++)
            {
                if (Primes5mod6_t[k-k_low_2])
                    ans.push_back((uint64_t)6*k-1);
                if (Primes1mod6_t[k-k_low_2])
                    ans.push_back((uint64_t)6*k+1);
            }
        }
    }
}
 
int main()
{
    // find vector ans of primes in range [6*k_start, 6*k_end+1]
    std::vector<uint64_t> ans;    

    int64_t  k_start=(int64_t)exp2(60)/6;          //k_start=0;
    int64_t  k_end=k_start+(int64_t)exp2(14)/6;    //k_end=(int64_t)1000000000/6;

    // dim step for first  segmented interval [0,sqrt(n)]
    int64_t  dim_step_1=L1_CACHE_SIZE-1;

    // dim step for second  segmented interval [sqrt(n),n]
    int64_t  dim_step_2=std::min(k_end-k_start+1,(int64_t)exp2(21)-1);

    if (k_end>k_start && k_start>=0)
        DoubleSegmentedSieve(k_start,k_end,dim_step_1,dim_step_2,ans);
  
    for (uint64_t p : ans)
    {
        std::cout << p << " ";
    }
    std::cout << std::endl;
    std::cout << ans.size()<< std::endl;
    
    return 0;
}