Prime sieve for large numbers

Question

I state that I am not an expert and that certainly the code can be improved.

I made this prime number sieve to be able to handle large numbers.

The basic idea is to write a number p=r+bW*k where r is the remainder of p divided by bW and bW the modulus. In practice it is a wheel sieve where multiples of the prime numbers divisors of bW are not stored. In this way, if n is the limit, we work with numbers smaller than n/bW and therefore we can increase bW to manage n ever larger ones.

bW can be increased by increasing n_PB

Of course, if for example n=10⁸ and bW=30, the sieve is designed to use n/bW=3333333 as a limit, so it counts the prime numbers up to 99999991.

The sieve is fast for example the execution time is 0.09s for n_PB=3 and input (0, 3333334) and less than 20s for n_PB=4 and input (0, 47619048) corresponds to n=10^10 and bW=210

///     This is a implementation of the bit wheel segmented sieve 

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

const int64_t PrimesBase[8]={2,3,5,7,11,13,17,19};
const int64_t n_PB = 3; // 3<= n_PB <=8

int64_t bW=1;
int64_t nR=0;
std::vector<int64_t> RW;
std::vector<int64_t> C_t;

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
};

int64_t Euclidean_Diophantine( int64_t coeff_a, int64_t  coeff_b)
{
    // return y in  Diophantine equation  coeff_a x + coeff_b y  = 1
    int64_t k=1;
    std::vector<int64_t> div_t; 
    std::vector<int64_t> rem_t;
    std::vector<int64_t> coeff_t;
    div_t.push_back(coeff_a);
    rem_t.push_back(coeff_b);
    coeff_t.push_back((int64_t)0);
    div_t.push_back((int64_t)div_t[0]/rem_t[0]);
    rem_t.push_back((int64_t)div_t[0]%rem_t[0]);
    coeff_t.push_back((int64_t)0);
    while (rem_t[k]>1)
    {
        k=k+1;
        div_t.push_back((int64_t)rem_t[k-2]/rem_t[k-1]);
        rem_t.push_back((int64_t)rem_t[k-2]%rem_t[k-1]);
        coeff_t.push_back((int64_t)0);
    }
    k=k-1;
    coeff_t[k]=-div_t[k+1];
    if (k>0)
        coeff_t[k-1]=(int64_t)1;
    while (k > 1)
    {
        k=k-1;
        coeff_t[k-1]=coeff_t[k+1];
        coeff_t[k]+=(int64_t)(coeff_t[k+1]*(-div_t[k+1]));
    }
    if (k==1)
        return (int64_t)(coeff_t[k-1]+coeff_t[k]*(-div_t[k]));
    else
        return (int64_t)(coeff_t[0]);
}

void get_wheel_constant(void)
{
    //get bW base wheel equal to bW=p1*p2*...*pn   with n=n_PB
    for(int64_t i=0; i<n_PB; i++)
        bW*=PrimesBase[i];
    //find reduct residue set
    std::vector<char> Remainder_t(bW,true); 
    for (int64_t i=0; i< n_PB; i++)
        for (int64_t j=PrimesBase[i];j< bW;j+=PrimesBase[i])
            Remainder_t[j]=false;

    for (int64_t j=2; j< bW; j++)
        if (Remainder_t[j]==true)
            RW.push_back(-bW+j);
    RW.push_back(1);
    nR=RW.size(); //nR=phi(bW)
    
    for (int64_t j=0; j<nR-2; j++)
        C_t.push_back(Euclidean_Diophantine(bW,-RW[j]));
    C_t.push_back(-1);
    C_t.push_back(1);   
}

int64_t get_mmin( int64_t k, int64_t ir, int64_t  jc)
{
    int64_t mmin=1;
    int64_t rW_t;
    if (ir==nR-1)
        mmin=0;
    if (jc==nR-1)
    {
        rW_t=RW[ir];
        mmin=0;                                 
    }
    else if (jc==ir)
    {
        rW_t=1;
        mmin=0;
    }
    else if (jc==nR-2)
    {
        rW_t=-bW-RW[ir];
        if (ir==nR-1)
            rW_t=-1;
    }
    else
    {
        rW_t=(C_t[jc]*(-RW[ir]))%bW;
        if(rW_t>1)
            rW_t-=bW;
    }
    mmin+=bW*k*k + k*(rW_t+RW[jc]) + (rW_t*RW[jc])/bW;
    return mmin;        
}

void segmented_bit_sieve_wheel(int64_t k_start,int64_t k_end)
{
    //count primes  from 1+bW*k_start to  1+bW*(k_end-1)
    int64_t count_p=0;

    int64_t segment_size_min=8191;

    get_wheel_constant();
    
    if (k_start<=PrimesBase[n_PB-1]/bW){
        count_p=n_PB;
        k_start=0;
    }
    if (k_end>1+PrimesBase[n_PB-1]/bW && k_end>bW && k_end>k_start && bW>=30){

        int64_t k_sqrt = (int64_t) std::sqrt(k_end/bW)+1;

        int64_t  nB=nR/8;
        int64_t segment_size=1; 
        int64_t p_mask_i=3;
        for (int64_t i=0; i<p_mask_i;i++)
            segment_size*=(bW+RW[i]);  
        while (segment_size<std::max(k_sqrt,segment_size_min) && p_mask_i<std::min(nR,(int64_t)7))
        {
            segment_size*=(bW+RW[p_mask_i]);  
            p_mask_i++;
        }
        
        int64_t segment_size_b=nB*segment_size;
        std::vector<uint8_t> Primes(nB+segment_size_b, 0xff);
        std::vector<uint8_t> Segment_i(nB+segment_size_b, 0xff);
        
        int64_t  pb,mb,ib,i,jb,j,k,kb;
        int64_t kmax = (int64_t) std::sqrt(segment_size/bW)+2;
        for (k =1; k  <= kmax; k++)
        {
            kb=k*nB;            
            for (jb = 0; jb<nB; jb++)
            {
                for (j = 0; j<8; j++)
                {
                    if(Primes[kb+jb] & (1 << j))
                    {
                        for (ib = 0; ib<nB; ib++)
                        {
                            for (i = 0; i<8; i++)
                            {
                                pb=nB*(bW*k+RW[j+jb*8]);                                
                                mb=nB*get_mmin( k, i+ib*8, j+jb*8);
                                for (; mb <= segment_size_b && mb>=0; mb +=pb )
                                    Primes[mb+ib] &= del_bit[i];
                                if (pb<nB*(bW+RW[p_mask_i]) && k_end>segment_size)
                                {
                                    mb-=segment_size_b;
                                    while (mb<0)
                                        mb+=pb;
                                    for (; mb <= segment_size_b; mb +=pb )
                                        Segment_i[mb+ib] &= del_bit[i];
                                }
                            }
                        }
                    }
                }
            }
        }
        if (k_start<segment_size)
        {
            for (kb = nB+nB*k_start; kb < std::min (nB+segment_size_b,nB*k_end); kb++)
                count_p+=bit_count[Primes[kb]]; 
        }
 
        if (k_end>segment_size) 
        {
            int64_t k_low,kb_low; 
            std::vector<uint8_t> Segment_t(nB+segment_size_b);
            k_low =segment_size;
            if (k_start>segment_size)
                k_low =(k_start/segment_size)*segment_size;
            
            for (; k_low < k_end; k_low += segment_size)
            {
                kb_low=k_low*nB;
                for (kb = 0; kb <(nB+segment_size_b); kb++)
                    Segment_t[kb]=Segment_i[kb]; 
                kmax=(std::min(segment_size,(int64_t)std::sqrt((k_low+segment_size)/bW)+2));
                j=p_mask_i;
                for(k=1; k<=kmax;k++)
                {
                    kb=k*nB;
                    for (jb = 0; jb<nB; jb++)
                    {
                        for (; j < 8; j++)
                        {
                            if (Primes[kb+jb]& (1 << j))
                            {
                                for (ib = 0; ib<nB; ib++)
                                {
                                    for (i = 0; i < 8; i++)
                                    {
                                        pb=bW*k+RW[j+jb*8];
                                        mb=-k_low+get_mmin(k, i+ib*8, j+jb*8);
                                        if (mb<0)
                                            mb=(mb%pb+pb)%pb;
                                        mb*=nB;
                                        pb*=nB;
                                        for (; mb <= segment_size_b; mb += pb)
                                            Segment_t[mb+ib] &= del_bit[i];
                                    }
                                }
                            }
                        }
                        j=0;
                    }
                }
                kb =nB+kb_low;
                if (k_start>k_low)
                    kb =nB+nB*k_start;
                for ( ; kb <std::min (kb_low+segment_size_b+nB,nB*k_end); kb++)
                    count_p+=bit_count[Segment_t[kb-kb_low]];
            }

        }
    }

    std::cout << " primes count= "  << count_p<< std::endl;
}

int main()
{
    //segmented_bit_sieve_wheel(k_start,k_end)
    //bW=PrimesBase[0]*PrimesBase[1]*...*PrimesBase[n_PB-1]
    // if n_PB=3 bW=30 if n_PB=4 bW=210 ... bW=2310 ... if n_PB=8 bW=9699690
    //k_start=n_start/bW and k_end=n_end/bW+1 find primes from 1+k_start*bW to 1+bW*(k_end-1)
    segmented_bit_sieve_wheel(0,3333334); 

    return 0;
}

To make it faster you have to use multithreading but I have no idea how to do it.

I thought that you can create m different blocks for multithreading in the second part of the code after for(; k_low < k_end; k_low += segment_size), using m different Segment_t vectors and increasing for each block k_low by m*segment_size.

It's a good idea?

Answer 1

Prefer C++ headers. <stdint.h> is a (deprecated) C header, that doesn't respect namespaces. Prefer the C++ version <cstdint>, which defines its identifiers in the std namespace.

I don't see a need for exactly-64-bit signed integers. All the arithmetic looks unsigned, so something like std::uint_fast64_t might be a better (and at least technically more portable) choice.

The variable names convey very little meaning to me - is there no way it can be clearer? The global variables in particular are a concern - not only do they need better names, but we should strive to avoid globals as much as possible, as they represent invisible coupling between functions - prefer to pass objects so that we can control that coupling.

We shouldn't be printing from within a computation function - separate that from the generation of primes. In particular, it's pointless benchmarking functions that do output.

Your speculation that dividing the work between threads is likely accurate, and I recommend you try that and compare the performance to the single-threaded version.