mersenneforum.org How to efficiently sieve for prime k-tuples
 Register FAQ Search Today's Posts Mark Forums Read

 2023-01-24, 15:52 #1 makis     "Makis Christou" Jun 2020 Cyprus B16 Posts 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 page that explains how rieMiner's algorithm works. I can say that I understand all the steps until optimizations come into play . This is my understanding so far: 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}$$ 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!
 2023-01-24, 18:32 #2 henryzz Just call me Henry     "David" Sep 2007 Liverpool (GMT/BST) 136448 Posts 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$$.
2023-01-25, 12:34   #3
makis

"Makis Christou"
Jun 2020
Cyprus

138 Posts

Quote:
 Originally Posted by henryzz 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$$.
Thanks for the explanation! It really helped me clear up my confusion. I think I got the basic idea working in a simple Python script. 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 .

2023-01-25, 16:29   #4
henryzz
Just call me Henry

"David"
Sep 2007
Liverpool (GMT/BST)

136448 Posts

Quote:
 Originally Posted by makis Thanks for the explanation! It really helped me clear up my confusion. I think I got the basic idea working in a simple Python script. 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 .
You might be interested in this thread: https://mersenneforum.org/showthread.php?t=16705 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.

 2023-01-25, 17:41 #5 User140242   Jul 2022 10011002 Posts I don't know if it's a good solution but this sieve wheel 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
 2023-01-28, 16:21 #6 User140242   Jul 2022 22×19 Posts 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 #include #include #include #include #include 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 RW; std::vector 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 div_t; std::vector rem_t; std::vector 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 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; j1) 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 &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 Primes(nB+segment_size_b, 0xff); std::vector 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=0; mb +=pb ) Primes[mb+ib] &= del_bit[i]; if (pbsegment_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_startsegment_size) { int64_t k_low,kb_low; std::vector 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; jbk_low) kb =nB+nB*k_start; for ( ; kb Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30}; search_k_tuples(0, 47619, Constellation); return 0; }
 2023-01-28, 18:52 #7 User140242   Jul 2022 22×19 Posts 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 #include #include #include #include #include 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 RW; std::vector 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 div_t; std::vector rem_t; std::vector 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 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; j1) 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 &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 Primes(nB+segment_size_b, 0xff); std::vector 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=0; mb +=pb ) Primes[mb+ib] &= del_bit[i]; if (pbsegment_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_startsegment_size) { int64_t k_low,kb_low; std::vector 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; jbk_low) k =1+k_start; for ( ; k Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30}; search_k_tuples(0, 47619, Constellation); return 0; } The sieve was not designed to do this, feedback is appreciated. Last fiddled with by User140242 on 2023-01-28 at 18:59
2023-01-29, 07:24   #8
makis

"Makis Christou"
Jun 2020
Cyprus

11 Posts

Quote:
 Originally Posted by User140242 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 #include #include #include #include #include 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 RW; std::vector 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 div_t; std::vector rem_t; std::vector 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 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; j1) 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 &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 Primes(nB+segment_size_b, 0xff); std::vector 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=0; mb +=pb ) Primes[mb+ib] &= del_bit[i]; if (pbsegment_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_startsegment_size) { int64_t k_low,kb_low; std::vector 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; jbk_low) k =1+k_start; for ( ; k Constellation = {0, 4, 10, 12, 18, 22, 24, 28, 30}; search_k_tuples(0, 47619, Constellation); return 0; }` The sieve was not designed to do this, feedback is appreciated.
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.

2023-01-29, 12:44   #9
User140242

Jul 2022

22·19 Posts

Quote:
 Originally Posted by makis 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.

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 i_c 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 bW = pm# is set, all multiples of pi 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 RW[i_c] + bW * (k_start + 1) is equivalent to T' + pm# * f + o and the value of i_c is searched finding the constellation between consecutive remainders given that the multiples of primes pi are excluded at the outset.

Last fiddled with by User140242 on 2023-01-29 at 13:40

 2023-03-01, 07:48 #10 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 2·5·7·13 Posts The code in Perl/ntheory in C (code) and C+GMP (code) 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. Last fiddled with by danaj on 2023-03-01 at 07:49

 Similar Threads Thread Thread Starter Forum Replies Last Post jzakiya jzakiya 111 2022-09-25 04:10 AddieJ Miscellaneous Math 15 2021-05-11 14:57 MattcAnderson Math 0 2020-06-16 17:45 Puzzle-Peter Software 156 2019-06-03 20:19 fivemack Math 27 2015-12-12 18:42

All times are UTC. The time now is 10:18.

Fri Mar 24 10:18:30 UTC 2023 up 218 days, 7:47, 0 users, load averages: 0.75, 0.97, 0.97