Created
August 7, 2022 13:51
-
-
Save user140242/efd065b221b02c548dce395a6ccd0c2f to your computer and use it in GitHub Desktop.
segmented sieve wheel base 30
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/// This is a implementation of the wheel segmented sieve with the basis {2,3,5}. | |
/// This algorithm uses (sqrt(n) *8/30) space. | |
#include <iostream> | |
#include <cmath> | |
#include <algorithm> | |
#include <vector> | |
#include <cstdlib> | |
#include <stdint.h> | |
const int64_t L1_CACHE_SIZE = 8192; | |
const std::vector<int64_t> C1={-22, -32, -28, -32, -28, -8, -8, -22, | |
-30, -18, -30, -30, -12, -30, -12, -18, | |
-34, -26, -16, -14, -34, -26, -14, -16, | |
-42, -42, -18, -12, -18, -18, -18, -12, | |
-46, -20, -34, -26, -10, -14, -20, -10, | |
-24, -36, -36, -24, -24, -6, -24, -6, | |
-36, -30, -24, -36, -30, -24, 0, 0, | |
-40, -38, -40, -20, -22, -20, -2, 2}; | |
const std::vector<int64_t> C2={0, 9, 7, 9, 7, 1, 1, 0, | |
6, 0, 8, 8, 1, 6, 1, 0, | |
9, 5, 0, 1, 9, 5, 1, 0, | |
15, 15, 1, 0, 3, 3, 1, 0, | |
18, 1, 10, 6, 0, 2, 1, 0, | |
1, 11, 11, 5, 5, 0, 1, 0, | |
10, 7, 4, 10, 7, 4, 0, 0, | |
13, 12, 13, 3, 4, 3, 0, 0}; | |
const std::vector<int64_t> RW={-23, -19, -17, -13, -11,-7 , -1, 1}; | |
const int64_t bW=30; | |
const int64_t nR=8; | |
void Sieve_Wheel_30(int64_t dim_step, std::vector<char> &Primes) | |
{ | |
int64_t m,mmin,i,j; | |
int64_t kmax=(int64_t)std::sqrt(dim_step/bW)+2; | |
for (int64_t k = 1; k <= kmax; k++) | |
{ | |
for (j = 0; j < nR; j++) | |
{ | |
if(Primes[k+j*(dim_step+1)]) | |
{ | |
for (i = 0; i < nR; i++) | |
{ | |
// Cancelling out all the multiples of bW*k+RW[j] | |
mmin=bW*k*k + k*C1[i*nR+j] + C2[i*nR+j]; | |
for ( m =mmin; m <= dim_step; m += bW*k+RW[j]) | |
Primes[m+i*(dim_step+1)] = false; | |
} | |
} | |
} | |
} | |
} | |
void Segmented_Sieve_Wheel_30(int64_t n) | |
{ | |
// counting primes <n n>bW | |
int64_t count_p=(n < 3) ? 0 : ((n == 3) ? 1 : ((n <5) ? 2 : 3)); | |
if (n>6) | |
{ | |
int64_t k_end=(n < bW) ? 2 :(int64_t)n/bW+1; | |
int64_t k_sqrt=(int64_t)std::sqrt(k_end/bW)+1; | |
int64_t dim_step,k; | |
int64_t k_low_2=0; | |
dim_step=std::max(L1_CACHE_SIZE-2,k_sqrt); | |
std::vector<char> Primes(nR*(dim_step + 1), true); | |
Sieve_Wheel_30(dim_step, Primes); | |
for (k = 1; k <= std::min(dim_step,k_end-1); k++) | |
for (int64_t i = 0; i < nR; i++) | |
if(Primes[k+i*(dim_step+1)]) | |
count_p++; | |
if (k==k_end && k<=dim_step && k>0) | |
for (int64_t i = 0; i<nR; i++) | |
if(Primes[k+i*(dim_step+1)] and RW[i]<(int64_t)(n%bW-bW)) | |
count_p++; | |
if (k_end>dim_step) | |
{ | |
int64_t i,j,kmax,m,mmin,p_j; | |
std::vector<char> Primes_t(nR*(dim_step + 1)); | |
for(k_low_2=dim_step;k_low_2<k_end;k_low_2+=dim_step) | |
{ | |
std::fill(Primes_t.begin(), Primes_t.end(), true); | |
kmax=std::min((int64_t)std::sqrt((k_low_2+dim_step)/bW)+2,dim_step); | |
for(k=1; k<=kmax;k++) | |
{ | |
for (j = 0; j < nR; j++) | |
{ | |
if (Primes[k+j*(dim_step+1)]) | |
{ | |
for (i = 0; i < nR; i++) | |
{ | |
p_j=bW*k+RW[j]; | |
mmin=-k_low_2+bW*k*k + k*C1[i*nR+j] + C2[i*nR+j]; | |
mmin=(mmin>=0)?mmin:(mmin%p_j+p_j)%p_j; | |
for (m =mmin; m <= dim_step; m += p_j) | |
Primes_t[m+i*(dim_step+1)] = false; | |
} | |
} | |
} | |
} | |
for ( k =1+k_low_2; k <=std::min (k_low_2+dim_step,k_end-1); k++) | |
for (i = 0; i < nR; i++) | |
if (Primes_t[k-k_low_2+i*(dim_step+1)]) | |
count_p++; | |
} | |
if (k==k_end && k-k_low_2<=dim_step && k-k_low_2>0) | |
for (i = 0; i<nR; i++) | |
if(Primes_t[k-k_low_2+i*(dim_step+1)] && RW[i]<(int64_t)(n%bW-bW)) | |
count_p++; | |
} | |
} | |
std::cout << " primes < " << n << ": "<< count_p<< std::endl; | |
} | |
int main() | |
{ | |
Segmented_Sieve_Wheel_30(100000000); | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment