2
\$\begingroup\$

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

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

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

bW can be increased by increasing n_PB

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

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

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

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

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

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

const int64_t del_bit[8] =
{
  ~(1 << 0),~(1 << 1),~(1 << 2),~(1 << 3),
  ~(1 << 4),~(1 << 5),~(1 << 6),~(1 << 7)
};

const int64_t bit_count[256] =
{
  0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
  4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
};

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

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

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

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

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

    int64_t segment_size_min=8191;

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

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

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

        }
    }

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

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

    return 0;
}

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

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

It's a good idea?

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

1 Answer 1

1
\$\begingroup\$

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

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

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

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

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

\$\endgroup\$
2
  • \$\begingroup\$ I used global variables just to be more compact and more readable. In reality the only variables with meaning besides the global ones are k_start and k_end, the others are simple variables to scroll the vector. The sieve is the same as the traditional sieve instead of p * p we use (r+ bW * k) * (s+ bW * k), with s chosen appropriately so that (r * s)%bW corresponds to the row, or residue, we are considering. Elements are stored in an array of nR=phi(bW) rows and k_end columns (which are segmented). To store a single column 8 * nB bits are used, with nB=nR/8, and then (uint8_t)*nB is used. \$\endgroup\$
    – user140242
    Sep 21 at 12:08
  • \$\begingroup\$ If (r+ bW * k) * (s+ bW * k) <(r * s)%bW+k_low * bW then mb is negative for this the sign is needed. \$\endgroup\$
    – user140242
    Sep 21 at 12:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.