Prime sieve for large numbers

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;
segment_size*=(bW+RW[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];
{
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));
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?

• Micro-review - prefer C++ header <cstdint>, which declares std::int64_t etc. Do you need exactly 64 bits here, or would the "fast" or "least" variant be more appropriate? Commented Sep 21, 2022 at 9:57
• @TobySpeight I'm sorry but I don't understand exactly what you are asking me. I know that the code can be improved, I hope it's clear how it works but I need to know if you can create a multithreading version in order to speed it up further. Commented Sep 21, 2022 at 10:28

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.

• I used global variables just to be more compact and more readable. In reality the only variables with meaning besides the global ones are k_start and k_end, the others are simple variables to scroll the vector. The sieve is the same as the traditional sieve instead of p * p we use (r+ bW * k) * (s+ bW * k), with s chosen appropriately so that (r * s)%bW corresponds to the row, or residue, we are considering. Elements are stored in an array of nR=phi(bW) rows and k_end columns (which are segmented). To store a single column 8 * nB bits are used, with nB=nR/8, and then (uint8_t)*nB is used. Commented Sep 21, 2022 at 12:08
• If (r+ bW * k) * (s+ bW * k) <(r * s)%bW+k_low * bW then mb is negative for this the sign is needed. Commented Sep 21, 2022 at 12:11
• I made a multithreaded version using uint64_t but I don't know if there is an improvement using std::uint_fast64_t. Commented Jan 21 at 17:02
• It's quite likely that your platform has std::uint64_t and std::uint_fast64_t the same type. But on other platforms, std::uint64_t may be absent, and you'd like your code to work there too. Commented Jan 21 at 17:06