Remarkable results with psyco and sieve of Eratosthenes

Discussion in 'Python' started by Steve Bergman, Nov 29, 2006.

  1. Just wanted to report a delightful little surprise while experimenting
    with psyco.
    The program below performs astonoshingly well with psyco.

    It finds all the prime numbers < 10,000,000

    Processor is AMD64 4000+ running 32 bit.

    Non psyco'd python version takes 94 seconds.

    psyco'd version takes 9.6 seconds.

    But here is the kicker. The very same algorithm written up in C and
    compiled with gcc -O3, takes 4.5 seconds. Python is runng half as fast
    as optimized C in this test!

    Made my day, and I wanted to share my discovery.

    BTW, can this code be made any more efficient?


    ============

    #!/usr/bin/python -OO
    import math
    import sys
    import psyco

    psyco.full()

    def primes():
    primes=[3]
    for x in xrange(5,10000000,2):
    maxfact = int(math.sqrt(x))
    flag=True
    for y in primes:
    if y > maxfact:
    break
    if x%y == 0:
    flag=False
    break
    if flag == True:
    primes.append(x)
    primes()
    Steve Bergman, Nov 29, 2006
    #1
    1. Advertising

  2. Steve Bergman

    Will McGugan Guest


    > #!/usr/bin/python -OO
    > import math
    > import sys
    > import psyco
    >
    > psyco.full()
    >
    > def primes():
    > primes=[3]
    > for x in xrange(5,10000000,2):
    > maxfact = int(math.sqrt(x))
    > flag=True
    > for y in primes:
    > if y > maxfact:
    > break
    > if x%y == 0:
    > flag=False
    > break
    > if flag == True:
    > primes.append(x)
    > primes()
    >


    Some trivial optimizations. Give this a whirl.

    def primes():
    sqrt=math.sqrt
    primes=[3]
    for x in xrange(5,10000000,2):
    maxfact = int(sqrt(x))
    for y in primes:
    if y > maxfact:
    primes.append(x)
    break
    if not x%y:
    break
    return primes

    --
    blog: http://www.willmcgugan.com
    Will McGugan, Nov 29, 2006
    #2
    1. Advertising

  3. Steve Bergman

    Will McGugan Guest

    Steve Bergman wrote:
    > Just wanted to report a delightful little surprise while experimenting
    > with psyco.
    > The program below performs astonoshingly well with psyco.
    >
    > It finds all the prime numbers < 10,000,000


    Actualy, it doesn't. You forgot 1 and 2.


    Will McGugan
    --
    blog: http://www.willmcgugan.com
    Will McGugan, Nov 29, 2006
    #3
  4. Steve Bergman

    Beliavsky Guest

    Will McGugan wrote:
    > Steve Bergman wrote:
    > > Just wanted to report a delightful little surprise while experimenting
    > > with psyco.
    > > The program below performs astonoshingly well with psyco.
    > >
    > > It finds all the prime numbers < 10,000,000

    >
    > Actualy, it doesn't. You forgot 1 and 2.


    The number 1 is not generally considered to be a prime number -- see
    http://mathworld.wolfram.com/PrimeNumber.html .
    Beliavsky, Nov 29, 2006
    #4
  5. Steve Bergman

    Will McGugan Guest

    Will McGugan, Nov 29, 2006
    #5
  6. Steve Bergman

    Guest

    > BTW, can this code be made any more efficient?

    I'm not sure, but the following code takes around 6 seconds on my
    1.2Ghz iBook. How does it run on your machine?

    def smallPrimes(n):
    """Given an integer n, compute a list of the primes < n"""

    if n <= 2:
    return []

    sieve = range(3, n, 2)
    top = len(sieve)
    for si in sieve:
    if si:
    bottom = (si*si - 3)//2
    if bottom >= top:
    break
    sieve[bottom::si] = [0] * -((bottom-top)//si)

    return [2]+filter(None, sieve)

    smallPrimes(10**7)




    >
    > ============
    >
    > #!/usr/bin/python -OO
    > import math
    > import sys
    > import psyco
    >
    > psyco.full()
    >
    > def primes():
    > primes=[3]
    > for x in xrange(5,10000000,2):
    > maxfact = int(math.sqrt(x))
    > flag=True
    > for y in primes:
    > if y > maxfact:
    > break
    > if x%y == 0:
    > flag=False
    > break
    > if flag == True:
    > primes.append(x)
    > primes()
    , Nov 29, 2006
    #6
  7. Will McGugan wrote:

    > Some trivial optimizations. Give this a whirl.



    I retimed and got 9.7 average for 3 runs on my version.

    Yours got it down to 9.2.

    5% improvement. Not bad.

    (Inserting '2' at the beginning doesn't seem to impact performance
    much.;-) )

    BTW, strictly speaking, shouldn't I be adding something to the floating
    point sqrt result, before converting to int, to allow for rounding
    error? If it is supposed to be 367 but comes in at 366.99999999, don't
    I potentially classify a composite as a prime?

    How much needs to be added?
    Steve Bergman, Nov 29, 2006
    #7
  8. wrote:
    > > BTW, can this code be made any more efficient?

    >
    > I'm not sure, but the following code takes around 6 seconds on my
    > 1.2Ghz iBook. How does it run on your machine?
    >
    >


    Hmm. Come to think of it, my algorithm isn't the sieve.

    Anyway, this is indeed fast as long as you have enough memory to handle
    it for the range supplied.

    It comes in at 1.185 seconds average on this box.

    Come to think of it, there is a supposedly highly optimized version of
    the sieve in The Python Cookbook that I've never bothered to actually
    try out. Hmmm...
    Steve Bergman, Nov 29, 2006
    #8
  9. Steve Bergman

    Guest

    On Nov 29, 6:59 pm, "Steve Bergman" <> wrote:
    > wrote:
    > > > BTW, can this code be made any more efficient?

    >
    > > I'm not sure, but the following code takes around 6 seconds on my
    > > 1.2Ghz iBook. How does it run on your machine?

    >
    > Hmm. Come to think of it, my algorithm isn't the sieve.


    Right. I guess the point of the sieve is that you don't have to spend
    any time
    finding that a given odd integer is not divisible by a given prime;
    all the
    multiplies are done up front, so you save all the operations
    corresponding to
    the case when x % y != 0 in your code. Or something.

    > Anyway, this is indeed fast as long as you have enough memory to handle
    > it for the range supplied.


    The sieve can be segmented, so that the intermediate space requirement
    for
    computing the primes up to n is O(sqrt(n)). (Of course you'll still
    need
    O(n/log n) space to store the eventual list of primes.) Then there
    are all sorts
    of bells and whistles (not to mention wheels) that you can add to
    improve the
    running time, most of which would considerably complicate the code.

    The book by Crandall and Pomerance (Primes: A Computational
    Perspective)
    goes into plenty of detail on all of this.

    Mark Dickinson
    , Nov 30, 2006
    #9
  10. On Wed, 29 Nov 2006 15:35:39 -0800, Steve Bergman wrote:

    > BTW, strictly speaking, shouldn't I be adding something to the floating
    > point sqrt result, before converting to int, to allow for rounding
    > error?


    If you don't mind doing no more than one unnecessary test per candidate,
    you can just add one to maxfact to allow for that. Or use round()
    rather than int(). Or don't convert it at all, just say:

    maxfact = math.sqrt(x)

    and compare directly to that.


    > If it is supposed to be 367 but comes in at 366.99999999, don't
    > I potentially classify a composite as a prime?


    Do you fear the math.sqrt() function is buggy? If so, all bets are off :)


    > How much needs to be added?


    No more than 1, and even that might lead you to sometimes performing an
    unnecessary test.



    --
    Steven.
    Steven D'Aprano, Nov 30, 2006
    #10
  11. At Wednesday 29/11/2006 20:35, Steve Bergman wrote:

    >BTW, strictly speaking, shouldn't I be adding something to the floating
    >point sqrt result, before converting to int, to allow for rounding
    >error? If it is supposed to be 367 but comes in at 366.99999999, don't
    >I potentially classify a composite as a prime?


    You could avoid sqrt using divmod (which gets the % result too); stop
    when quotient<=divisor.
    But this approach creates a tuple and then unpacks it, so you should
    time it to see if there is an actual speed improvement.


    --
    Gabriel Genellina
    Softlab SRL

    __________________________________________________
    Correo Yahoo!
    Espacio para todos tus mensajes, antivirus y antispam ¡gratis!
    ¡Abrí tu cuenta ya! - http://correo.yahoo.com.ar
    Gabriel Genellina, Nov 30, 2006
    #11
  12. Will McGugan wrote:
    > > #!/usr/bin/python -OO
    > > import math
    > > import sys
    > > import psyco
    > >
    > > psyco.full()
    > >
    > > def primes():
    > > primes=[3]
    > > for x in xrange(5,10000000,2):
    > > maxfact = int(math.sqrt(x))
    > > flag=True
    > > for y in primes:
    > > if y > maxfact:
    > > break
    > > if x%y == 0:
    > > flag=False
    > > break
    > > if flag == True:
    > > primes.append(x)
    > > primes()
    > >

    >
    > Some trivial optimizations. Give this a whirl.
    >
    > def primes():
    > sqrt=math.sqrt
    > primes=[3]
    > for x in xrange(5,10000000,2):
    > maxfact = int(sqrt(x))
    > for y in primes:
    > if y > maxfact:
    > primes.append(x)
    > break
    > if not x%y:
    > break
    > return primes


    You can also save an attribute lookup for append; just add
    append = primes.append
    outside of the loop and replace primes.append(x) with append(x)
    That should cut down a few fractions of second.

    George
    George Sakkis, Nov 30, 2006
    #12
  13. Steve Bergman

    Guest

    George Sakkis:
    > You can also save an attribute lookup for append; just add
    > append = primes.append
    > outside of the loop and replace primes.append(x) with append(x)
    > That should cut down a few fractions of second.


    We were talking about Psyco, and I think with Psyco (just released for
    Py 2.5, BTW) such tricks are less useful.

    Bye,
    bearophile
    , Nov 30, 2006
    #13
  14. In article <>,
    Steve Bergman wrote:
    >BTW, can this code be made any more efficient?


    >def primes():
    > primes=[3]
    > for x in xrange(5,10000000,2):
    > maxfact = int(math.sqrt(x))
    > flag=True
    > for y in primes:
    > if y > maxfact:
    > break

    [...]

    You can omit the call to math.sqrt if you test this instead.

    y*y > x

    in place of if y > maxfact: .

    Pka
    Pekka Karjalainen, Nov 30, 2006
    #14
  15. Pekka Karjalainen wrote:

    > You can omit the call to math.sqrt if you test this instead.
    >
    > y*y > x
    >
    > in place of if y > maxfact: .


    Or use

    sqrt = lambda x: x ** .5

    Cheers,

    --
    Klaus Alexander Seistrup
    http://klaus.seistrup.dk/
    Klaus Alexander Seistrup, Nov 30, 2006
    #15
  16. Steve Bergman

    Klaas Guest

    Klaus Alexander Seistrup wrote:
    > Pekka Karjalainen wrote:
    >
    > > You can omit the call to math.sqrt if you test this instead.
    > >
    > > y*y > x
    > >
    > > in place of if y > maxfact: .

    >
    > Or use
    >
    > sqrt = lambda x: x ** .5


    Test it:

    $ python -m timeit -s "from math import sqrt" "sqrt(5.6)"
    1000000 loops, best of 3: 0.445 usec per loop
    $ python -m timeit -s "sqrt = lambda x: x**.5" "sqrt(5.6)"
    1000000 loops, best of 3: 0.782 usec per loop

    Note that this overhead is almost entirely in function calls; calling
    an empty lambda is more expensive than a c-level sqrt:

    $ python -m timeit -s "sqrt = lambda x: x" "sqrt(5.6)"
    1000000 loops, best of 3: 0.601 usec per loop

    Just math ops:
    $ python -m timeit -s "x = 5.6" "x*x"
    10000000 loops, best of 3: 0.215 usec per loop
    $ python -m timeit -s "x = 5.6" "x**.5"
    1000000 loops, best of 3: 0.438 usec per loop

    Of course, who knows that psyco does with this under the hood.

    -Mike
    Klaas, Nov 30, 2006
    #16
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. jzakiya
    Replies:
    3
    Views:
    539
    jzakiya
    Jul 19, 2008
  2. jzakiya
    Replies:
    4
    Views:
    465
    user923005
    Jun 13, 2008
  3. Daniel Baird

    golfing Eratosthenes

    Daniel Baird, Aug 2, 2006, in forum: Ruby
    Replies:
    13
    Views:
    227
    Chad Perrin
    Aug 2, 2006
  4. jzakiya
    Replies:
    6
    Views:
    178
    Michael Ulm
    Jun 11, 2008
  5. Sebastian

    Java 8 Streams and Eratosthenes

    Sebastian, Jun 4, 2013, in forum: Java
    Replies:
    22
    Views:
    638
    Sebastian
    Jun 18, 2013
Loading...

Share This Page