[Tutor] primes - sieve of odds

[C Smith]

On Wednesday, Mar 23, 2005, at 04:00 America/Chicago, 
tutor-request at python.org wrote:

> Now it would be fine to have an *equally fast*
> infinite prime number generator.
>
> Has anybody any suggestions?
>

I think when you add the condition of it being an infinite generator, 
you are changing the rules and can't end up with something as fast.  
What makes the sieve so fast is that we are generating all the primes 
in a given range using a very simple "strike out" method.  In the 
infinite generator scenario, you will get all the primes up to n with 
the sieve and can yield those back, but at some point you will have to 
stop and get "the next one." To get the next one you will have to pick 
a range in which a next prime is likely to be found. The previous 
primes can be used to clear out this range.

What follows is an attempt based on the previous tutor-evolved sieve 
that extends the range in which to find the next prime by a factor of 2 
over the last known prime. A similar algorithm is on ASPN (I bellieve), 
under

Space-efficient version of sieve of Eratosthenes.
D. Eppstein, May 2004

It's slower than the sieve but about twice as fast as Eppstein's up to 
1,000,000.  I'll leave it to someone else to see when it finally blows 
up :-) The output of primes was checked through 1,000,000 against 
Eppstein's generator without error.

/c

###
def esieve():
'''extended sieve generator that returns primes on each call. When the 
end of the
existing list is reached, more primes are sought in the range that 
extends to
twice the last known prime. The known primes are used to clear out this 
extended
range using the Sieve of Eratosthenes.'''

     # get started with primes less than 100; you need at least [2, 3]
     primes=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
         41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
     i=0
     while 1:
         yield primes[i]
         i+=1
         if i==len(primes): #extend range
             minn=primes[-1]+2 #this is an odd number
             maxx=2*minn+1 #there should be a prime in this range; +1 
makes it odd
             sqrt=int(maxx**.5) #don't use primes bigger than this
             length = (maxx-minn)//2 #this is used for computing the 
crossing out None values
             nums=range(minn,maxx+1,2) #here is the range in which a 
primes will be found
             for p in primes:
                 if p>sqrt:break
                 j=minn%p #find out where the striking out should start
                 if j<>0:
                     j=(p-j) #if evens were present, this is where to 
start, but...
                     if (minn+j)%2==0:j+=p #if it lands on an even, go 
to the next one
                     j//=2 #and now account for the fact that evens 
aren't present
                 nums[j::p]=[None]*(1+(length-j)//p) #cross out 
multiples of p
             x=filter(None,nums) #clean up the range
             assert len(x)>0 #there should be a prime in here, but check 
anyway
             primes.extend(x) #add these to the existing primes and 
continue yielding
###