mersenneforum.org  

Go Back   mersenneforum.org > Great Internet Mersenne Prime Search > Software

Reply
 
Thread Tools
Old 2023-01-24, 15:52   #1
makis
 
makis's Avatar
 
"Makis Christou"
Jun 2020
Cyprus

1110 Posts
Question 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:
  1. Pick Primorial Number \(m\), based on desired digit number (the bigger the better but I am sure practical limitations come into play)
  2. Calculate offset \(o\) (based on Primorial and desired pattern, this is pre-calculated in rieMiner)
  3. Calculate \(T^{'} = T + p_m \# - (T \; mod \; p_m \#)\)
  4. 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!
makis is offline   Reply With Quote
Old 2023-01-24, 18:32   #2
henryzz
Just call me Henry
 
henryzz's Avatar
 
"David"
Sep 2007
Liverpool (GMT/BST)

3×2,011 Posts
Default

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\).
henryzz is online now   Reply With Quote
Old 2023-01-25, 12:34   #3
makis
 
makis's Avatar
 
"Makis Christou"
Jun 2020
Cyprus

11 Posts
Smile

Quote:
Originally Posted by henryzz View Post
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 .
makis is offline   Reply With Quote
Old 2023-01-25, 16:29   #4
henryzz
Just call me Henry
 
henryzz's Avatar
 
"David"
Sep 2007
Liverpool (GMT/BST)

3·2,011 Posts
Default

Quote:
Originally Posted by makis View Post
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.
henryzz is online now   Reply With Quote
Old 2023-01-25, 17:41   #5
User140242
 
Jul 2022

22·19 Posts
Default

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

Last fiddled with by User140242 on 2023-01-25 at 18:21
User140242 is offline   Reply With Quote
Old 2023-01-28, 16:21   #6
User140242
 
Jul 2022

22·19 Posts
Default

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;
}
User140242 is offline   Reply With Quote
Old 2023-01-28, 18:52   #7
User140242
 
Jul 2022

22·19 Posts
Default

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;
}
The sieve was not designed to do this, feedback is appreciated.

Last fiddled with by User140242 on 2023-01-28 at 18:59
User140242 is offline   Reply With Quote
Old 2023-01-29, 07:24   #8
makis
 
makis's Avatar
 
"Makis Christou"
Jun 2020
Cyprus

138 Posts
Unhappy

Quote:
Originally Posted by User140242 View Post
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;
}
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.
makis is offline   Reply With Quote
Old 2023-01-29, 12:44   #9
User140242
 
Jul 2022

22·19 Posts
Default

Quote:
Originally Posted by makis View Post
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
User140242 is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
On the Infinity of Twin Prime and other K-tuples jzakiya jzakiya 111 2022-09-25 04:10
How to get to a prime from the previous prime number (sieve of Eratosthenes rediscovered) AddieJ Miscellaneous Math 15 2021-05-11 14:57
prime k-tuples MattcAnderson Math 0 2020-06-16 17:45
How do you efficiently sieve for prime 3/4-tuples? Puzzle-Peter Software 156 2019-06-03 20:19
Efficiently finding a linear progression in data fivemack Math 27 2015-12-12 18:42

All times are UTC. The time now is 14:30.


Tue Feb 7 14:30:21 UTC 2023 up 173 days, 11:58, 1 user, load averages: 0.57, 0.90, 1.09

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2023, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.

≠ ± ∓ ÷ × · − √ ‰ ⊗ ⊕ ⊖ ⊘ ⊙ ≤ ≥ ≦ ≧ ≨ ≩ ≺ ≻ ≼ ≽ ⊏ ⊐ ⊑ ⊒ ² ³ °
∠ ∟ ° ≅ ~ ‖ ⟂ ⫛
≡ ≜ ≈ ∝ ∞ ≪ ≫ ⌊⌋ ⌈⌉ ∘ ∏ ∐ ∑ ∧ ∨ ∩ ∪ ⨀ ⊕ ⊗ 𝖕 𝖖 𝖗 ⊲ ⊳
∅ ∖ ∁ ↦ ↣ ∩ ∪ ⊆ ⊂ ⊄ ⊊ ⊇ ⊃ ⊅ ⊋ ⊖ ∈ ∉ ∋ ∌ ℕ ℤ ℚ ℝ ℂ ℵ ℶ ℷ ℸ 𝓟
¬ ∨ ∧ ⊕ → ← ⇒ ⇐ ⇔ ∀ ∃ ∄ ∴ ∵ ⊤ ⊥ ⊢ ⊨ ⫤ ⊣ … ⋯ ⋮ ⋰ ⋱
∫ ∬ ∭ ∮ ∯ ∰ ∇ ∆ δ ∂ ℱ ℒ ℓ
𝛢𝛼 𝛣𝛽 𝛤𝛾 𝛥𝛿 𝛦𝜀𝜖 𝛧𝜁 𝛨𝜂 𝛩𝜃𝜗 𝛪𝜄 𝛫𝜅 𝛬𝜆 𝛭𝜇 𝛮𝜈 𝛯𝜉 𝛰𝜊 𝛱𝜋 𝛲𝜌 𝛴𝜎𝜍 𝛵𝜏 𝛶𝜐 𝛷𝜙𝜑 𝛸𝜒 𝛹𝜓 𝛺𝜔