wheel segmented sieve

wheel_segmented_sieve.cpp

///     This is a implementation of the  wheel segmented sieve for counting primes.

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

void Sieve_Wheel_Base_6(int64_t  dim_step, std::vector<char> &Primes5mod6, std::vector<char> &Primes1mod6)
{
    // Find 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 Wheel_Segmented_Sieve( int64_t max_base_wheel, int64_t  k_end,int64_t dim_step_2)
{
    int64_t  sqrt_max=(int64_t)std::sqrt(6*k_end+1);
    int64_t  dim_step=(int64_t)std::sqrt(sqrt_max)/6;
    int64_t  k_low,k_low_2;
    int64_t  num_step=0;
    int64_t  base_modulo=1;
    int64_t  len_primebase=0;
    std::vector<int64_t> Primes_i; 
    
    // Find all the primes >3 <= sqrt(sqrt_max)+1
    std::vector<char> Primes5mod6(dim_step + 1, true);
    std::vector<char> Primes1mod6(dim_step + 1, true);
    Sieve_Wheel_Base_6(dim_step, Primes5mod6, Primes1mod6);

    // Append in Primes_i  primes in range [0, sqrt(sqrt_max)+1]
    Primes_i.push_back(2);
    Primes_i.push_back(3);
    for (int64_t k = 1; k <= dim_step; k++)
    {
        if (Primes5mod6[k])
            Primes_i.push_back(6*k-1);
        if (Primes1mod6[k])
            Primes_i.push_back(6*k+1);
    } 
    
    int64_t  len_primes_dim_step=Primes_i.size();


    std::vector<int64_t> PrimesBase({2,3,5,7,11,13});
    
	//get modulo base wheel equal to p1*p2*...*pn <=max_base_wheel  with n=len_primebase
    while(len_primebase<PrimesBase.size()&&(base_modulo*PrimesBase[len_primebase]<=std::min(max_base_wheel,sqrt_max)))
    {
        base_modulo*=PrimesBase[len_primebase];
        len_primebase++;
    }
	
	//inizialize count=utilized number primes base 
    int64_t   countp=len_primebase;
	
	//find possible remainder of base module
    std::vector<char> Remainder_i_t(base_modulo+1,true); 
    for (int64_t i=0; i< len_primebase;i++)
        for (int64_t j=PrimesBase[i]*PrimesBase[i];j< base_modulo+1;j+=PrimesBase[i])
             Remainder_i_t[j]=false;
    std::vector<int64_t> Remainder_i;
    for (int64_t j=PrimesBase[len_primebase-1]+1;j< base_modulo+1;j++)
        if (Remainder_i_t[j]==true)
            Remainder_i.push_back(-base_modulo+j);
    Remainder_i.push_back(1);
    int64_t  len_Remainder_i=Remainder_i.size();

    //temporary vectors for coefficients calculating the initial mmin index 
	//equivalent to p * p in traditional sieve
    std::vector<int64_t> Remainder_low(len_Remainder_i);
    std::vector<int64_t> Remainder_low_div(len_Remainder_i);
    std::vector<int64_t> Remainder_low2(len_Remainder_i);
    std::vector<int64_t> Remainder_low_div2(len_Remainder_i);
    
    //dim_step_1 dimension
    int64_t  dim_step_1=(int64_t)sqrt_max/base_modulo;
    int64_t  mod_sqrt_max=(int64_t)sqrt_max%base_modulo;
    if (mod_sqrt_max>1)
    {
        mod_sqrt_max-=base_modulo;
        dim_step_1++;
    }

    // temp vectors find primes  remainder Remainder_i[i2] in a single step [low,low+dim_step_2] in the interval [sqrt_max+c+1,k_end]
    std::vector<char> Primes_f(dim_step_2 + 1, true);
   //temp vectors find primes  remainder Remainder_i[i]in range [0,sqrt_max]
    std::vector<char> Primes_t(dim_step_1 + 1,true);  

    int ck_ans=0;
    for(k_low_2=dim_step_1-1;k_low_2<=(int64_t) 2+k_end/(base_modulo/6);k_low_2+=dim_step_2)
    {
        int64_t  mmin,mmin1,p_k_i,k_i,p_j,k_j,m_j;
        int64_t i,remaindertemp;
        for (int64_t i2=0; i2< len_Remainder_i;i2++)
        {
			//coefficients calculating for the initial mmin1 index 
            for (int64_t j=0; j< len_Remainder_i;j++)
            {
                for (int64_t k=0; k< len_Remainder_i;k++)
                {
                    int64_t mod1=0;
                    remaindertemp=Remainder_i[k]*Remainder_i[j]%base_modulo;
                    if (remaindertemp>1)
                    {
                        remaindertemp-=base_modulo;
                        mod1=1;
                    }
                    if (remaindertemp==Remainder_i[i2])
                    {
                        Remainder_low2[j]=Remainder_i[k]+Remainder_i[j];
                        Remainder_low_div2[j]=Remainder_i[k]*Remainder_i[j]/base_modulo+mod1;
                        break;
                    }
                }
            }
            for (i=0; i< len_Remainder_i;i++)
            {
                //coefficients calculating for the initial mmin index 
				if (i==i2)
				{
					Remainder_low=Remainder_low2;
					Remainder_low_div=Remainder_low_div2;
				}
				else
				{
					for (int64_t j=0; j< len_Remainder_i;j++)
					{
						for (int64_t k=0; k< len_Remainder_i;k++)
						{
							int64_t mod1=0;
							remaindertemp=Remainder_i[k]*Remainder_i[j]%base_modulo;
							if (remaindertemp>1)
							{
								remaindertemp-=base_modulo;
								mod1=1;
							}
							if (remaindertemp==Remainder_i[i])
							{
								Remainder_low[j]=Remainder_i[k]+Remainder_i[j];
								Remainder_low_div[j]=Remainder_i[k]*Remainder_i[j]/base_modulo+mod1;
								break;
							}
						}
					}
				}
				// Cancelling out all the multiples of p_j = m_j + base_modulo *k_j in Primes_t
                for(int64_t j=len_primebase; j<len_primes_dim_step;j++)
                {
                    p_j=(int64_t)Primes_i[j];
                    k_j=(int64_t)p_j/base_modulo;
                    m_j=(int64_t)p_j%base_modulo;
                    if (m_j>1)
                    {
                         m_j-=base_modulo;
                         k_j++;
                    }
					int64_t i_m_j=0;
					for (; i_m_j< len_Remainder_i;i_m_j++)
					if (m_j==Remainder_i[i_m_j])
					    break;
					mmin=(int64_t)base_modulo*k_j*k_j+k_j*Remainder_low[i_m_j]+Remainder_low_div[i_m_j];
					if (mmin <=dim_step_1 && mmin>=0)    
					for (int64_t m =mmin; m <= dim_step_1; m += p_j)
						Primes_t[m] = false;
				}
				// Counting primes remainder Remainder_i[i]  of range [0,sqrt_max]
				if (ck_ans==0)
				{
					int64_t kans;
					for (kans=1; kans < dim_step_1; kans++)
						if (Primes_t[kans])
							countp++;
					if (Remainder_i[i]<=mod_sqrt_max && kans <= dim_step_1)
						if (Primes_t[kans])
							countp++;
				}
				// Cancelling out all the multiples of p_k_i=Remainder_i[i]+k_i*base_modulo in Primes_f
				int64_t  kimax=(int64_t)std::sqrt((k_low_2+dim_step_2)/base_modulo)+2;
				kimax=std::min(kimax,dim_step_1);
				for (int64_t k_i=1; k_i <= kimax; k_i++)
				{
					p_k_i=Remainder_i[i]+k_i*base_modulo;
					if (Primes_t[k_i])
					{
						mmin1=(int64_t)base_modulo*k_i*k_i+k_i*Remainder_low2[i]+Remainder_low_div2[i];
						mmin1=(mmin1>=k_low_2)?mmin1-k_low_2:((mmin1-k_low_2)%p_k_i+p_k_i)%p_k_i;
						if (mmin1>=0 && mmin1 <= dim_step_2)
							for (int64_t m1 =mmin1; m1 <= dim_step_2; m1 += p_k_i)
								Primes_f[m1] = false;
					}
				}
				std::fill(Primes_t.begin(), Primes_t.end(), true);
			}
			ck_ans=1;
			// Counting primes remainder Remainder_i[i2] in step dim_step_2 of  [sqrt_max,6k_end+1]
			int64_t kans2=1+k_low_2;
			if (kans2==dim_step_1)
			{
				if (Remainder_i[i]>mod_sqrt_max)
					if (Primes_f[1])
						countp++;
				kans2++;
			}
			for (; kans2 <=std::min(k_low_2+dim_step_2, (int64_t)k_end/(base_modulo/6)) ; kans2++)
                if (Primes_f[kans2-k_low_2])
                    countp++;
			if (kans2==(int64_t)k_end/(base_modulo/6)+1)
			{
				for (;kans2 <=std::min(k_low_2+dim_step_2, (int64_t)k_end/(base_modulo/6)+2); kans2++)
					if ((int64_t)(kans2*base_modulo/6+Remainder_i[i2]/6)<=k_end+1)
						if (Primes_f[kans2-k_low_2])
							countp++;
			}
			std::fill(Primes_f.begin(), Primes_f.end(), true);
		}
    }
    std::cout << countp << " ";
    std::cout << std::endl;
}
 
int main()
{
    int64_t  k_end=(int64_t)100000000/6;
    
    int64_t  max_base_wheel=30;  //max value sqrt(6*k_end+1) modulo base wheel=p1*p2*...*pn <=max_base_wheel

    // dimension of step for second  segmented interval [sqrt(n),n]
    int64_t  dim_step_2=(int64_t)exp2(21)-1;

    if ((int64_t)std::sqrt(std::sqrt(6*k_end+1))/6>=(int64_t)16)
        Wheel_Segmented_Sieve(max_base_wheel,k_end,dim_step_2);
   
    return 0;
}