Rollback to Revision 1

Source Link

edited Mar 23, 2020 at 8:00

13.5k
12
54
119

#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,ih,il;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;
            }
            tl=i*30+res[j];
            tl=tl*tl;
            th=tl/30;th=i; // high
            tl=tl%30;tl=res[j]; // low
            // number30*th+res[tl] is 30*th+res[tl]composite
            while(th<t1){
                for(k=0;k<8;++k){th+=i;
                    if(tl==res[k]){tl+=res[j];
                if(tl>=30){
        primes[th]|=1<<k; // not a prime
        tl-=30;
                break;
    th+=1;
                }
        // adding prime to self
    }
            if(th>=t){
    th+=2*i;
                tl+=2*res[j];break;
                if(tl>=60){
 } // exceeds bound
                tl-=60;for(k=0;k<8;++k){
                    th+=2;if(tl==res[k]){
                }else if(tl>=30){
            primes[th]|=1<<k; // not a prime
    tl-=30;
                    th+=1;break;
                } // adding number to self}
                }
            }
        }
    }
    // 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;
}

Multiples of 3 without optimization: 7.69
Multiples of 30 without optimization: 1428.5542
Multiples of 3 with optimization: 4.00
Multiples of 30 with optimization: 47.9432

Edit: Added 1201ProgramAlarm's comment on starting on i**2 and it is still slower, however the program is still slower, perhaps the increments will have to change to avoid checking if the number mod 30 is 1,7,11,13...

#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,ih,il;
    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<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;
            }
            tl=i*30+res[j];
            tl=tl*tl;
            th=tl/30; //high
            tl=tl%30; // low
            // number is 30*th+res[tl]
            while(th<t){
                for(k=0;k<8;++k){
                    if(tl==res[k]){
                        primes[th]|=1<<k; // not a prime
                        break;
                    }
                }
                th+=2*i;
                tl+=2*res[j];
                if(tl>=60){
                    tl-=60;
                    th+=2;
                }else if(tl>=30){
                    tl-=30;
                    th+=1;
                } // adding number to self
                
            }
        }
    }
    // 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;
}

Multiples of 3 without optimization: 7.69
Multiples of 30 without optimization: 14.55
Multiples of 3 with optimization: 4.00
Multiples of 30 with optimization: 4.94

Edit: Added 1201ProgramAlarm's comment on starting on i**2 and it is still slower, however the program is still slower, perhaps the increments will have to change to avoid checking if the number mod 30 is 1,7,11,13...

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

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

updated code based on comment

Source Link

edited Mar 23, 2020 at 4:36

Ariana

153
7

#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;tl,ih,il;
    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;tl=i*30+res[j];
 /           tl=tl*tl;
            th=tl/30; //high
            tl=res[j];tl=tl%30; // low
            // 30*th+res[tl]number is composite30*th+res[tl]
            while(1th<t){
                th+=i;for(k=0;k<8;++k){
                tl+=res[j];
    if(tl==res[k]){
            if(tl>=30){
            primes[th]|=1<<k; // not a prime
    tl-=30;
                    th+=1;break;
                } // adding prime to self}
                if(th>=t){}
                th+=2*i;
    break;
            tl+=2*res[j];
    } // exceeds bound         if(tl>=60){
                for(k=0;k<8;++k){    tl-=60;
                    if(tl==res[k]){th+=2;
                }else if(tl>=30){
            primes[th]|=1<<k; // not a prime
    tl-=30;
                    break;th+=1;
                } // adding number }to self
                }
            }
        }
    }
    // 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;
}

Multiples of 3 without optimization: 7.69
Multiples of 30 without optimization: 2814.4255
Multiples of 3 with optimization: 4.00
Multiples of 30 with optimization: 74.3294

Edit: Added 1201ProgramAlarm's comment on starting on i**2 and it is still slower, however the program is still slower, perhaps the increments will have to change to avoid checking if the number mod 30 is 1,7,11,13...

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

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

#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,ih,il;
    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<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;
            }
            tl=i*30+res[j];
            tl=tl*tl;
            th=tl/30; //high
            tl=tl%30; // low
            // number is 30*th+res[tl]
            while(th<t){
                for(k=0;k<8;++k){
                    if(tl==res[k]){
                        primes[th]|=1<<k; // not a prime
                        break;
                    }
                }
                th+=2*i;
                tl+=2*res[j];
                if(tl>=60){
                    tl-=60;
                    th+=2;
                }else if(tl>=30){
                    tl-=30;
                    th+=1;
                } // adding number to self
                
            }
        }
    }
    // 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;
}

Multiples of 3 without optimization: 7.69
Multiples of 30 without optimization: 14.55
Multiples of 3 with optimization: 4.00
Multiples of 30 with optimization: 4.94

Edit: Added 1201ProgramAlarm's comment on starting on i**2 and it is still slower, however the program is still slower, perhaps the increments will have to change to avoid checking if the number mod 30 is 1,7,11,13...

Source Link

asked Mar 22, 2020 at 16:15

Ariana

153
7

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

performance c sieve-of-eratosthenes

Stack Exchange Network

Return to Question

Sieve of Eratosthenes optimisation