[Tutor] primes (generator)

Gregor Lingl glingl at aon.at
Sat Mar 19 01:04:57 CET 2005


Hi all of you!

Many thanks for your remarks, ideas and experiments.

When I posted my input yesterday:

[x for x in range(2,100) if not [y for y in range(2,x) if x%y==0]]

my intention was indeed to find a "shortest" a program, disregarding
efficiency. Thus range instead of xrange etc.

(About one year ago I was involved in a "coding contest" at Linux-Tag
in Karlsruhe, where the goal was to achive a certain output with a
shortest program ;-) Strange. I used Python, but I didn't win :-(  )

Clearly there are several criteria concerning the quality of programs
as e. g. conciseness, clarity, efficiency, beauty(!) etc. which may not
be achieved optimally all together ...

I definitely intend to sum up your contibutions (for me) later, but for
now (because of lack of time) I'll only present one more idea,
induced by Danny's (unfortunately buggy - because of recursion overflow) 
contribution:

Thought it would be nice to use Python's new generator-expressions
instead of itertool module. First idea:

##################################################################
# infinite prime generator
# induced by ideas of Danny Yoo,
# which themselves were induced by material of SICP

import math

def prime(n):
     return [x for x in range(2,1+int(math.sqrt(n))) if n%x==0]==[]

def infrange(n):
     while True:
	yield n
	n+=1

primes = (p for p in infrange(2) if prime(p))

###############  The generator primes works like this:

 >>> primes.next()
2
 >>> primes.next()
3
 >>> primes.next()
5
 >>> [primes.next() for i in range(10)]
[7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
 >>> [primes.next() for i in range(100)]
[43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 
113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 
193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 
359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 
443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 
541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617]
 >>>
 >>> ### You may also use it this way:
 >>> for p in primes:
	if p > 1000: break
	print p,

	
619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 
743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 
863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997


With a fresh generator

[primes.next() for p in xrange(10000)]

produces (on my machine) a list of the first 10000 primes (upto 104729)
within 4.68 secs.

Now for a bit of optimization: instead of building a list of
all divisors use a loop and get out as soon as possible:

So I replace prime by:

def prime(n):
     for d in range(2,1+int(math.sqrt(n))):
         if n%d == 0: return  False
     return True

This makes the generator nearly 4 times as fast. I get my list in
1.33 secs.

Lets optimize a bit more: test only 2 and odd numbers:

def infrange(n):
     yield 2
     n=3
     while True:
	yield n
	n+=2

and do a similar thing for searching for divisors:

def prime(n):
     if n==2: return True
     for d in [2]+ range(3,1+int(math.sqrt(n)),2):
         if n%d == 0: return  False
     return True

Now the first 10000  primes get computed in 0.67 secs.
Again double speed. Together seven times as fast as
first version. But IMO at the cost of clarity
and conciseness.

Maybe you have amendments (of course a matter of taste)
to this code. I'd be interested in them ..

Maybe sombody likes to compare these algorithms to the
ones mentioned before in this thread. I would be
interested in the results.

Regards,
Gregor



Danny Yoo schrieb:
> 
> On Thu, 17 Mar 2005, Gregor Lingl wrote:
> 
> 
>>Hi!
>>Who knows a more concise or equally concise but more efficient
>>expression, which returns the same result as
...
> ######
> 
>>>>from itertools import ifilter, count
>>>>
>>>>def notDivisibleBy(n):
> 
> ...     def f(x):
> ...         return x % n != 0
> ...     return f
> ...
> 
>>>>def sieve(iterable):
> 
> ...     nextPrime = iterable.next()
> ...     yield nextPrime
> ...     restOfPrimes = sieve(ifilter(notDivisibleBy(nextPrime), iterable))
> ...     for prime in restOfPrimes:
> ...         yield prime
> ...
> 
>>>>primes = sieve(count(2))
>>>>primes.next()
> 
> 2
> 
>>>>primes.next()
...
>>>>primes.next()
> 
> 23
> ######
> 
> The neat thing about this code is that it produces an infinite iterator of
> prime numbers.  The original code in Scheme is itself quite concise and
> quite nice.  It is described here:
> 
> http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-24.html#call_footnote_Temp_455
> 
> 
> Best of wishes to you!
> 
> 
> 

-- 
Gregor Lingl
Reisnerstrasse 3/19
A-1030 Wien

Telefon: +43 1 713 33 98
Mobil:   +43 664 140 35 27

Autor von "Python für Kids"
Website: python4kids.net


More information about the Tutor mailing list