How to efficiently sieve for prime k-tuples
 

So I have been trying to understand how efficient sieving for prime k-tuples works for quite a while but now actually dedicated some time to understand the maths behind it. I am looking at this [URL="https://riecoin.dev/en/Mining_Algorithm"]page[/URL] that explains how rieMiner's algorithm works. I can say that I understand all the steps until optimizations come into play :smile:. This is my understanding so far:
[LIST=1][*]Pick Primorial Number [$]m[/$], based on desired digit number (the bigger the better but I am sure practical limitations come into play)[*]Calculate offset [$]o[/$] (based on Primorial and desired pattern, this is pre-calculated in rieMiner)[*]Calculate [$]T^{'} = T + p_m \# - (T \; mod \; p_m \#)[/$][*]Check if [$]T_f = T^{'} + p_m\# \cdot f + o[/$] is a base prime of a constellation of the desired pattern for different values of f up to an arbitrary [$]f_{max}[/$][/LIST]
This is the basic algorithm without any optimizations. The idea behind the optimized version is to efficiently eliminate some f's without requiring the expensive primality test. Based on the explanation in the "Basis of the current implementation" section, the idea to me at least, is that after picking a primorial number [$]P[/$] (which has changed from [$]m[/$] for some reason?), and calculating the offset [$]o[/$], we sieve out the factors [$]T_{f_p}[/$] that have an f of the form [$]f_p = ((p-((T^{'} + o) \; mod \; p)) {p_P\#}^{-1}) \; mod \; p[/$] as well as their multiples i.e. [$]T_{p \cdot i + f_p}[/$].

My confusion comes from the fact that I am not sure how we start sieving the different [$]f_p[/$] or what the p's represent in this case. Are they primes in our prime table? Are they something else? Also is the inverse of the primorial we use to compute [$]f_p[/$] different for each [$]f_p[/$]? Or is it the initial one we pick? I am not really a number theory guy so if anyone is kind enough to provide some pseudo-code that would be very helpful. Thanks!

Like you I am not finding their notation easy to follow so I will start from scratch.

The aim is to sieve the polynomial [$]k*n\#+o+c_i[/$] where o is an offset and c_i are the prime constellation terms(maybe 0,2,6 etc) for a range of k.
To sieve each prime p it is necessary to work out the first k where [$]k*n\#+o+c_i[/$] is divisible by p for each [$]c_i[/$] or [$]k*n\#+o+c_i \equiv 0\text{ mod }p[/$].
To do this we solve for k, [$]k \equiv (-o-c_i)/n\#\text{ mod }p[/$]. The division is done by calculating the modular inverse of n# mod p.
Once we have this first k we can find more by adding p until k>max k.

In their notation, each [$]f_p[/$] is a prime that they are sieving(p above). The modular inverse is different for each [$]f_p[/$].

[QUOTE=henryzz;623393]Like you I am not finding their notation easy to follow so I will start from scratch.

The aim is to sieve the polynomial [$]k*n\#+o+c_i[/$] where o is an offset and c_i are the prime constellation terms(maybe 0,2,6 etc) for a range of k.
To sieve each prime p it is necessary to work out the first k where [$]k*n\#+o+c_i[/$] is divisible by p for each [$]c_i[/$] or [$]k*n\#+o+c_i \equiv 0\text{ mod }p[/$].
To do this we solve for k, [$]k \equiv (-o-c_i)/n\#\text{ mod }p[/$]. The division is done by calculating the modular inverse of n# mod p.
Once we have this first k we can find more by adding p until k>max k.

In their notation, each [$]f_p[/$] is a prime that they are sieving(p above). The modular inverse is different for each [$]f_p[/$].[/QUOTE]

Thanks for the explanation! It really helped me clear up my confusion. I think I got the basic idea working in a simple [URL="https://pastebin.com/1i9HnyRF"]Python script[/URL]. Currently my sieving implementation takes too long since I am using a dictionary and not an efficiently implemented boolean array. Regardless, I am getting the same number of prime k-tuples with a fraction of the primality tests :banana:.

[QUOTE=makis;623443]Thanks for the explanation! It really helped me clear up my confusion. I think I got the basic idea working in a simple [URL="https://pastebin.com/1i9HnyRF"]Python script[/URL]. Currently my sieving implementation takes too long since I am using a dictionary and not an efficiently implemented boolean array. Regardless, I am getting the same number of prime k-tuples with a fraction of the primality tests :banana:.[/QUOTE]

You might be interested in this thread: [url]https://mersenneforum.org/showthread.php?t=16705[/url] particularly later in the thread with polysieve. Polysieve can be persuaded to sieve this form by inputting the primorial, offset and constellation manually.
Recently I have been working on a primorial variation which simplifies this a bit. The original polysieve should perform just as well as my new version, although I have been experimenting with adding buckets to the second stage sieve(I believe the majority of time was in cache misses). More testing is needed before I release it.
I imagine rieMiner is likely to perform better as it considers SIMD, which polysieve doesn't.

I don't know if it's a good solution but this [URL="https://codereview.stackexchange.com/questions/279803/prime-sieve-for-large-numbers"]sieve wheel [/URL]uses modular arithmetic so p = r + bW * k with k>0

bW is the size of the wheel and r is the possible remainder

example if bW=210 then r is one of these values:

RW[48]={-199, -197, -193, -191, -187, -181, -179, -173, -169, -167, -163, -157, -151, -149, -143, -139, -137, -131, -127, -121, -113, -109, -107, -103, -101, -97, -89, -83, -79, -73, -71, -67, -61, -59, -53, -47, -43, -41, -37, -31, -29, -23, -19, -17, -13, -11, -1, 1};

In the sieve it is possible to increase the size of the wheel to search for larger numbers and the search is done by setting kmin and kmmax

for bW=210 for each value of k six bytes are used in order to have 48 bits corresponding to each value of the remainder.

If a delta vector is defined

delta_R[48]={0, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2};

once the constellation has been identified, it is possible to verify if the corresponding bits are set to 1.

The sieve in the link was created to count the prime numbers but it can be easily modified and once the index of the first term of the constellation has been identified

for example indicated with ic the index of the first term of the constellation and k_tuples is the number of elements, then it is possible to use in the second part this code to see if all the bits of the constellation are set to 1

[CODE]

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_tuples=0;
for (ib = ic/8; ib < nB; ib++)
for (i = ic%8; i < 8 && i+8*ib < ic + k_tuples; i++)
if(Segment_t[kb - kb_low + ib]& (1 << i))
count_tuples++;
if (count_tuples==k_tuples)
for (ib = ic/8; ib < nB; ib++)
for (i = ic%8; i < 8 && i+8*ib < ic + k_tuples; i++)
std::cout<<" , "<<RW[i+ib*8]+(kb+k_low)*bW;
}
[/CODE]

It needs to be perfected but this is a simple example:

[CODE]
/// This is a implementation of the bit wheel segmented sieve for search k-tuples

#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 = 4; // 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 search_k_tuples(int64_t k_start,int64_t k_end, std::vector<int> &Constellation)
{

int64_t count_tuple=0;
int64_t k_tuple = Constellation.size();

int64_t segment_size_min=8191;

get_wheel_constant();

if (k_start<=PrimesBase[n_PB-1]/bW)
k_start=0;

//search Constellation
int64_t i_c = nR;
for (int64_t i = 0; i < nR - k_tuple; i++)
{
int64_t j = 1;
while ((int)(RW[i + j] - RW[i]) == Constellation[j] && i + j < nR - 1 && j < k_tuple -1)
j++;
if (j == k_tuple - 1 && (int)(RW[i+j] - RW[i]) == Constellation[j])
{
i_c = i;
break;
}
}

if (i_c < nR)
{
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_tuple=0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Primes[kb + ib]& (1 << i))
count_tuple++;
i = 0 ;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + (kb/nB) * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

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_tuple = 0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Segment_t[kb - kb_low + ib] & (1 << i))
count_tuple++;
i = 0;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c/8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + (kb/nB + k_low) * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

}
}
}
else
std::cout << "invalid constellation" << std::endl;
}

int main()
{
std::vector<int> Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30};
search_k_tuples(0, 47619, Constellation);

return 0;
}
[/CODE]

I can't edit anymore Sorry but I noticed an error

[CODE]
/// This is a implementation of the bit wheel segmented sieve for search k-tuples

#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 = 4; // 3<= n_PB <=8

//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

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

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 search_k_tuples(int64_t k_start,int64_t k_end, std::vector<int> &Constellation)
{

int64_t count_tuple=0;
int64_t k_tuple = Constellation.size();

int64_t segment_size_min=8191;

get_wheel_constant();

if (k_start<=PrimesBase[n_PB-1]/bW)
k_start=0;

//search Constellation
int64_t i_c = nR;
for (int64_t i = 0; i < nR - k_tuple; i++)
{
int64_t j = 1;
while ((int)(RW[i + j] - RW[i]) == Constellation[j] && i + j < nR - 1 && j < k_tuple -1)
j++;
if (j == k_tuple - 1 && (int)(RW[i+j] - RW[i]) == Constellation[j])
{
i_c = i;
break;
}
}

if (i_c < nR)
{
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 (k = 1+k_start; k < std::min (1+segment_size,k_end); k++)
{
count_tuple=0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Primes[k*nB + ib]& (1 << i))
count_tuple++;
i = 0 ;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + k * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

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;
}
}
k =1+k_low;
if (k_start>k_low)
k =1+k_start;
for ( ; k <std::min (k_low+segment_size+1,k_end); k++)
{
count_tuple = 0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Segment_t[k*nB - kb_low + ib] & (1 << i))
count_tuple++;
i = 0;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c/8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + k * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

}
}
}
else
std::cout << "invalid constellation" << std::endl;
}

int main()
{
std::vector<int> Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30};
search_k_tuples(0, 47619, Constellation);

return 0;
}
[/CODE]

The sieve was not designed to do this, feedback is appreciated.

[QUOTE=User140242;623667]I can't edit anymore Sorry but I noticed an error

[CODE]
/// This is a implementation of the bit wheel segmented sieve for search k-tuples

#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 = 4; // 3<= n_PB <=8

//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

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

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 search_k_tuples(int64_t k_start,int64_t k_end, std::vector<int> &Constellation)
{

int64_t count_tuple=0;
int64_t k_tuple = Constellation.size();

int64_t segment_size_min=8191;

get_wheel_constant();

if (k_start<=PrimesBase[n_PB-1]/bW)
k_start=0;

//search Constellation
int64_t i_c = nR;
for (int64_t i = 0; i < nR - k_tuple; i++)
{
int64_t j = 1;
while ((int)(RW[i + j] - RW[i]) == Constellation[j] && i + j < nR - 1 && j < k_tuple -1)
j++;
if (j == k_tuple - 1 && (int)(RW[i+j] - RW[i]) == Constellation[j])
{
i_c = i;
break;
}
}

if (i_c < nR)
{
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 (k = 1+k_start; k < std::min (1+segment_size,k_end); k++)
{
count_tuple=0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Primes[k*nB + ib]& (1 << i))
count_tuple++;
i = 0 ;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + k * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

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;
}
}
k =1+k_low;
if (k_start>k_low)
k =1+k_start;
for ( ; k <std::min (k_low+segment_size+1,k_end); k++)
{
count_tuple = 0;
i = i_c % 8;
for (ib = i_c / 8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
if(Segment_t[k*nB - kb_low + ib] & (1 << i))
count_tuple++;
i = 0;
}
if (count_tuple == k_tuple)
{
i = i_c % 8;
for (ib = i_c/8; ib < nB; ib++)
{
for (; i < 8 && i + 8 * ib < i_c + k_tuple; i++)
std::cout << " - " << RW[i + ib * 8] + k * bW;
i = 0;
}
std::cout << std::endl;
}
}
}

}
}
}
else
std::cout << "invalid constellation" << std::endl;
}

int main()
{
std::vector<int> Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30};
search_k_tuples(0, 47619, Constellation);

return 0;
}
[/CODE]

The sieve was not designed to do this, feedback is appreciated.[/QUOTE]

Thanks for modifying your sieve for prime k-tuples but I am not sure I can understand your code, at least easily. Is this showing something different than the idea explained above by henryzz? I am firstly trying to understand the theory before jumping to (optimized) implementations.

[QUOTE=makis;623718]Thanks for modifying your sieve for prime k-tuples but I am not sure I can understand your code, at least easily. Is this showing something different than the idea explained above by henryzz? I am firstly trying to understand the theory before jumping to (optimized) implementations.[/QUOTE]


I tried to understand the algorithm in the link you posted and in fact the operation is slightly different but it should be equivalent.

In the example I posted only one[B] i_c[/B] index is searched but there could be more than one you only have to make a small change using a for loop.

To understand how it works, once[B] bW = p[SUB]m[/SUB]#[/B] is set, all multiples of p[SUB]i[/SUB] are automatically excluded and the constellation is searched among the possible remainders which are in number phi(bW) where phi is Euler's totient function.

In practice [B]RW[i_c] + bW * (k_start + 1)[/B] is equivalent to [B]T' + p[SUB]m[/SUB]# * f + o[/B] and the value of [B]i_c[/B] is searched finding the constellation between [B] consecutive remainders[/B] given that the multiples of primes p[SUB]i[/SUB] are excluded at the outset.

The code in Perl/ntheory in C ([URL="https://github.com/danaj/Math-Prime-Util/blob/master/sieve_cluster.c"]code[/URL]) and C+GMP ([URL="https://github.com/danaj/Math-Prime-Util-GMP/blob/2594cbd9a4d56643746bd24418b81ae86e5bea8e/gmp_main.c#L2084"]code[/URL]) does fairly efficient sieving for prime k-tuples. There are some simple Perl example scripts including one that does parallel searching.

A simple test was hundreds of times faster than the C++ sieve above, but that's not really fair given it does primality tests after filtering. Which of course for any reasonable size you will have to do, as you can't expect to do full sieves of 25+ digit numbers.

I post with the thought that perhaps the OP or others can find something of use.