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