Improving pure Python prime sieve by recurrence formula Improving pure Python prime sieve by recurrence formula python python

Improving pure Python prime sieve by recurrence formula


You could do a wheel optimisation. Multiples of 2 and 3 aren't primes, so dont store them at all. You could then start from 5 and skip multiples of 2 and 3 by incrementing in steps of 2,4,2,4,2,4 etc..

Below is a C++ code for it. Hope this helps.

void sieve23(){    int lim=sqrt(MAX);    for(int i=5,bit1=0;i<=lim;i+=(bit1?4:2),bit1^=1)    {        if(!isComp[i/3])        {            for(int j=i,bit2=1;;)            {                j+=(bit2?4*i:2*i);                bit2=!bit2;                if(j>=MAX)break;                isComp[j/3]=1;            }        }    }}


If you may decide you are going to C++ to improve the speed, I ported the Python sieve to C++. The full discussion can be found here: Porting optimized Sieve of Eratosthenes from Python to C++.

On Intel Q6600, Ubuntu 10.10, compiled with g++ -O3 and with N=100000000 this takes 415 ms.

#include <vector>#include <boost/dynamic_bitset.hpp>// http://vault.embedded.com/98/9802fe2.htm - integer square rootunsigned short isqrt(unsigned long a) {    unsigned long rem = 0;    unsigned long root = 0;    for (short i = 0; i < 16; i++) {        root <<= 1;        rem = ((rem << 2) + (a >> 30));        a <<= 2;        root++;        if (root <= rem) {            rem -= root;            root++;        } else root--;    }    return static_cast<unsigned short> (root >> 1);}// https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188// https://stackoverflow.com/questions/5293238/porting-optimized-sieve-of-eratosthenes-from-python-to-c/5293492template <class T>void primesbelow(T N, std::vector<T> &primes) {    T i, j, k, sievemax, sievemaxroot;    sievemax = N/3;    if ((N % 6) == 2) sievemax++;    sievemaxroot = isqrt(N)/3;    boost::dynamic_bitset<> sieve(sievemax);    sieve.set();    sieve[0] = 0;    for (i = 0; i <= sievemaxroot; i++) {        if (sieve[i]) {            k = (3*i + 1) | 1;            for (j = k*k/3; j < sievemax; j += 2*k) sieve[j] = 0;            for (j = (k*k+4*k-2*k*(i&1))/3; j < sievemax; j += 2*k) sieve[j] = 0;        }    }    primes.push_back(2);    primes.push_back(3);    for (i = 0; i < sievemax; i++) {        if (sieve[i]) primes.push_back((3*i+1)|1);    }}