# Sieve of Eratosthenes optimisation

I'm implementing Sieve of Eratosthenes by working with multiples of 30 and comparing it to multiples of 3 from a previous answer

code for multiples of 30:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(){
const unsigned int res[8] = {1,7,11,13,17,19,23,29};
const unsigned int N = 1000000000;
unsigned int i,j,k,th,tl;
u_int8_t *primes = calloc(N/30+1,sizeof(char));
// 0 is taken to be prime while 1 composite(opposite from the code for multiples of 3)
//jth bit of primes[i]: 30*i+res[j]
primes[0] = '\x01'; // initialize with 1 is not prime and the others are prime
unsigned int ub = sqrt(N)/30+1;
unsigned int t = N/30+1;
for(i=0;i<ub;++i){
for(j=0;j<8;++j){
//current number is i*30+res[j]
if(primes[i]>>j&1){// jth bit is set to 1
continue;
}
th=i; // high
tl=res[j]; // low
// 30*th+res[tl] is composite
while(1){
th+=i;
tl+=res[j];
if(tl>=30){
tl-=30;
th+=1;
} // adding prime to self
if(th>=t){
break;
} // exceeds bound
for(k=0;k<8;++k){
if(tl==res[k]){
primes[th]|=1<<k; // not a prime
break;
}
}
}
}
}
// counting primes
k=3; // 2,3,5
for(i=0;i<t-1;++i){
for(j=0;j<8;++j){
if(primes[i]>>j&1){
continue;
}
++k;
}
}
for(j=0;j<8;++j){
if(primes[i]>>j&1){
continue;
}
if(i*30+res[j]>N){
break;
}
++k;
}
printf("Number of primes equal or less than %d: %d\n",N,k);
free(primes);
return 0;
}


Timing both variants locally(with -O3 and without compiler optimization), this variant seems to perform worse than the one using multiples of 3:

Multiples of 3 without optimization: 7.69
Multiples of 30 without optimization: 28.42
Multiples of 3 with optimization: 4.00
Multiples of 30 with optimization: 7.32


looking at the output of -O3 for both programs, the compiler only unrolls the loop and hardcodes some computation(i.e. sqrt(N)) and that's basically it, so either taking multiples of 30 is slower theoretically or the implementation is slower, which is more likely to be the case.

Is there any way that this code can be optimized or is some better way to go about writing the sieve for multiples of 30?

--code for multiples of 3 used as comparison--

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

int main(void){
unsigned int N = 1000000000;
unsigned int arraySize = (N/24 + 1);
uint32_t *primes = malloc(arraySize);

// The bits in primes follow this pattern:
//
// Bit 0 = 5, bit 1 = 7, bit 2 = 11, bit 3 = 13, bit 4 = 17, etc.
//
// For even bits, bit n represents 5 + 6*n
// For odd  bits, bit n represents 1 + 6*n
memset(primes , 0xff, arraySize);

int sqrt_N = sqrt(N);
for(int i = 5; i <= sqrt_N; i += 4) {
int iBitNumber = i / 3 - 1;
int iIndex = iBitNumber >> 5;
int iBit   = 1 << (iBitNumber & 31);
if ((primes[iIndex] & iBit) != 0) {
int increment = i+i;
for (int j = i * i; j < N; j += increment) {
int jBitNumber = j / 3 - 1;
int jIndex = jBitNumber >> 5;
int jBit   = 1 << (jBitNumber & 31);

primes[jIndex] &= ~jBit;

j += increment;
if (j >= N)
break;

jBitNumber = j / 3 - 1;
jIndex = jBitNumber >> 5;
jBit   = 1 << (jBitNumber & 31);

primes[jIndex] &= ~jBit;

// Skip multiple of 3.
j += increment;
}
}
i += 2;
iBit <<= 1;
if ((primes[iIndex] & iBit) != 0) {
int increment = i+i;
for (int j = i * i; j < N; j += increment) {
int jBitNumber = j / 3 - 1;
int jIndex = jBitNumber >> 5;
int jBit   = 1 << (jBitNumber & 31);

primes[jIndex] &= ~jBit;

// Skip multiple of 3.
j += increment;

j += increment;
if (j >= N)
break;

jBitNumber = j / 3 - 1;
jIndex = jBitNumber >> 5;
jBit   = 1 << (jBitNumber & 31);

primes[jIndex] &= ~jBit;
}
}
}

// Initial count includes 2, 3.
int count=2;
for (int i=5;i<N;i+=6) {
int iBitNumber = i / 3 - 1;
int iIndex = iBitNumber >> 5;
int iBit   = 1 << (iBitNumber & 31);
if (primes[iIndex] & iBit) {
count++;
}
iBit <<= 1;
if (primes[iIndex] & iBit) {
count++;
}
}
printf("%d\n", count);

free(primes);
return 0;
}

• (I don't see multiples of 3 - should that be 6?) – greybeard Mar 22 '20 at 22:00
• @greybeard yea it kinda is in multiples of 6 (6n-1 or 6n+1) but was written as multiples of 3 in the prev ans so kinda carried it here – Ariana Mar 22 '20 at 23:31
• Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. – Mast Mar 23 '20 at 8:00

When the original multiples-of-3 code finds a prime, it starts setting bits with the square of that value (for (int j = i * i). Your multiples-of-30 code does not do this, and can waste a lot of time marking numbers "not prime" that have already been so marked. As the new prime gets larger, this will consume a growing amount of time.

Using 1201ProgramAlarm's comment, the code ran slightly faster but was still slower than taking multiples of 3(effectively 6). However finding if the current composite number is in the multplicative group mod 30 seemed to be the one that took up a long time and could have been memoisation. This is done by precomputing how many times one adds the prime to reach the next composite that needs to be marked and also precomputing the squares(since the counting starts from i*i).

Improved code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(){
const unsigned int res[8] = {1,7,11,13,17,19,23,29};
const unsigned int N = 1000000000;
unsigned int nextelem[8][8][2]={}; //res[i]+nextelem[i][j][0]*res[j]=nextelem[i][j][1]
unsigned int startk[8]={}; //res[i]*2=res[startk[j]]
unsigned int i,j,k,ii,jj,kk;
u_int8_t *primes = calloc(N/30+1,sizeof(char));
//jth bit of primes[i]: 30*i+res[j]
primes[0] = '\x01';
unsigned int ub = sqrt(N)/30+1;
unsigned int t = N/30+1;
for(i=0;i<8;++i){// genning nextelem
for(j=0;j<8;++j){
for(k=2;k<30;k+=2){
for(ii=0;ii<8;++ii){
if(res[ii]==(res[i]+k*res[j])%30){
break;
}
}
if(ii!=8){
nextelem[i][j][0]=k;
nextelem[i][j][1]=ii;
break;
}
}
}
}
for(i=0;i<8;++i){// genning startk
for(j=0,k=(res[i]*res[i])%30;j<8;++j){
if(res[j]==k){
startk[i]=j;
break;
}
}
}
for(i=0;i<ub;++i){
for(j=0;j<8;++j){
//current number is i*30+res[j]
if(primes[i]>>j&1){// jth bit is set to 1
continue;
}
// we start from the square and go up, have a lookup table to figure how much to increment
ii=i*30+res[j];
jj=ii*ii;
k=startk[j];
while(jj<N){
primes[jj/30]|=1<<k; // jj not a prime
jj+=nextelem[k][j][0]*ii;
k=nextelem[k][j][1];
}
}
}
// counting primes
k=3; // 2,3,5
for(i=0;i<t-1;++i){
for(j=0;j<8;++j){
if(primes[i]>>j&1){
continue;
}
++k;
}
}
for(j=0;j<8;++j){
if(primes[i]>>j&1){
continue;
}
if(i*30+res[j]>N){
break;
}
++k;
}
printf("Number of primes equal or less than %d: %d\n",N,k);
free(primes);
return 0;
}


Just a nice property the multiplicative group is actually isomorphic to Z2xZ4=(11)x(7) so startk only consists of either res[0]=1 or res[5]=7^2=19

Timings with -O3:

Multiples of 30: 2.70
Multiples of 3: 3.57


which is a nice improvement.