Ultimate Prime Sieve -- Sieve Of Zakiya (SoZ)

Discussion in 'Python' started by jzakiya, Jun 13, 2008.

  1. jzakiya

    jzakiya Guest

    This is to announce the release of my paper "Ultimate Prime Sieve --
    Sieve of Zakiiya (SoZ)" in which I show and explain the development of
    a class of Number Theory Sieves to generate prime numbers. I used
    Ruby 1.9.0-1 as my development environment on a P4 2.8 Ghz laptop.

    You can get the pdf of my paper and Ruby and Python source from here:

    http://www.4shared.com/dir/7467736/97bd7b71/sharing.html

    Below is a sample of one of the simple prime generators. I did a
    Python version of this in my paper (see Python source too). The Ruby
    version below is the minimum array size version, while the Python has
    array of size N (I made no attempt to optimize its implementation,
    it's to show the method).

    class Integer
    def primesP3a
    # all prime candidates > 3 are of form 6*k+1 and 6*k+5
    # initialize sieve array with only these candidate values
    # where sieve contains the odd integers representatives
    # convert integers to array indices/vals by i = (n-3)>>1 =
    (n>>1)-1
    n1, n2 = -1, 1; lndx= (self-1) >>1; sieve = []
    while n2 < lndx
    n1 +=3; n2 += 3; sieve[n1] = n1; sieve[n2] = n2
    end
    #now initialize sieve array with (odd) primes < 6, resize array
    sieve[0] =0; sieve[1]=1; sieve=sieve[0..lndx-1]

    5.step(Math.sqrt(self).to_i, 2) do |i|
    next unless sieve[(i>>1) - 1]
    # p5= 5*i, k = 6*i, p7 = 7*i
    # p1 = (5*i-3)>>1; p2 = (7*i-3)>>1; k = (6*i)>>1
    i6 = 6*i; p1 = (i6-i-3)>>1; p2 = (i6+i-3)>>1; k = i6>>1
    while p1 < lndx
    sieve[p1] = nil; sieve[p2] = nil; p1 += k; p2 += k
    end
    end
    return [2] if self < 3
    [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
    end
    end

    def primesP3(val):
    # all prime candidates > 3 are of form 6*k+(1,5)
    # initialize sieve array with only these candidate values
    n1, n2 = 1, 5
    sieve = [False]*(val+6)
    while n2 < val:
    n1 += 6; n2 += 6; sieve[n1] = n1; sieve[n2] = n2
    # now load sieve with seed primes 3 < pi < 6, in this case just 5
    sieve[5] = 5

    for i in range( 5, int(ceil(sqrt(val))), 2) :
    if not sieve: continue
    # p1= 5*i, k = 6*i, p2 = 7*i,
    p1 = 5*i; k = p1+i; p2 = k+i
    while p2 <= val:
    sieve[p1] = False; sieve[p2] = False; p1 += k; p2 += k
    if p1 <= val: sieve[p1] = False

    primes = [2,3]
    if val < 3 : return [2]
    primes.extend( i for i in range(5, val+(val&1), 2) if sieve )

    return primes

    Now to generate an array of the primes up to some N just do:

    Ruby: 10000001.primesP3a
    Python: primesP3a(10000001)

    The paper presents benchmarks with Ruby 1.9.0-1 (YARV). I would love
    to see my various prime generators benchmarked with optimized
    implementations in other languages. I'm hoping Python gurus will do
    better than I, though the methodology is very very simple, since all I
    do is additions, multiplications, and array reads/writes.

    Have fun with the code. ;-)

    Jabari Zakiya
    jzakiya, Jun 13, 2008
    #1
    1. Advertising

  2. jzakiya

    jzakiya Guest

    On Jun 13, 1:12 pm, jzakiya <> wrote:
    > This is to announce the release of my paper "Ultimate Prime Sieve --
    > Sieve of Zakiiya (SoZ)" in which I show and explain the development of
    > a class of Number Theory Sieves to generate prime numbers.   I used
    > Ruby 1.9.0-1 as my development environment on a P4 2.8 Ghz laptop.
    >
    > You can get the pdf of my paper and Ruby and Python source from here:
    >
    > http://www.4shared.com/dir/7467736/97bd7b71/sharing.html
    >
    > Below is a sample of one of the simple prime generators. I did a
    > Python version of this in my paper (see Python source too).  The Ruby
    > version below is the minimum array size version, while the Python has
    > array of size N (I made no attempt to optimize its implementation,
    > it's to show the method).
    >
    > class Integer
    >    def primesP3a
    >       # all prime candidates > 3 are of form  6*k+1 and 6*k+5
    >       # initialize sieve array with only these candidate values
    >       # where sieve contains the odd integers representatives
    >       # convert integers to array indices/vals by  i = (n-3)>>1 =
    > (n>>1)-1
    >       n1, n2 = -1, 1;  lndx= (self-1) >>1;  sieve = []
    >       while n2 < lndx
    >          n1 +=3;   n2 += 3;   sieve[n1] = n1;  sieve[n2] = n2
    >       end
    >       #now initialize sieve array with (odd) primes < 6, resize array
    >       sieve[0] =0;  sieve[1]=1;  sieve=sieve[0..lndx-1]
    >
    >       5.step(Math.sqrt(self).to_i, 2) do |i|
    >          next unless sieve[(i>>1) - 1]
    >          # p5= 5*i,  k = 6*i,  p7 = 7*i
    >          # p1 = (5*i-3)>>1;  p2 = (7*i-3)>>1;  k = (6*i)>>1
    >          i6 = 6*i;  p1 = (i6-i-3)>>1;  p2 = (i6+i-3)>>1;  k = i6>>1
    >          while p1 < lndx
    >              sieve[p1] = nil;  sieve[p2] = nil;  p1 += k;  p2 += k
    >          end
    >       end
    >       return [2] if self < 3
    >       [2]+([nil]+sieve).compact!.map {|i| (i<<1) +3 }
    >    end
    > end
    >
    > def primesP3(val):
    >     # all prime candidates > 3 are of form  6*k+(1,5)
    >     # initialize sieve array with only these candidate values
    >     n1, n2 = 1, 5
    >     sieve = [False]*(val+6)
    >     while  n2 < val:
    >         n1 += 6;   n2 += 6;  sieve[n1] = n1;   sieve[n2] = n2
    >     # now load sieve with seed primes 3 < pi < 6, in this case just 5
    >     sieve[5] = 5
    >
    >     for i in range( 5, int(ceil(sqrt(val))), 2) :
    >        if not sieve:  continue
    >        #  p1= 5*i,  k = 6*i,  p2 = 7*i,
    >        p1 = 5*i;  k = p1+i;  p2 = k+i
    >        while p2 <= val:
    >           sieve[p1] = False;  sieve[p2] = False;  p1 += k;  p2 += k
    >        if p1 <= val:  sieve[p1] = False
    >
    >     primes = [2,3]
    >     if val < 3 : return [2]
    >     primes.extend( i for i in range(5, val+(val&1), 2)  if sieve )
    >
    >     return primes
    >
    > Now to generate an array of the primes up to some N just do:
    >
    > Ruby:    10000001.primesP3a
    > Python: primesP3a(10000001)
    >
    > The paper presents benchmarks with Ruby 1.9.0-1 (YARV).  I would love
    > to see my various prime generators benchmarked with optimized
    > implementations in other languages.  I'm hoping Python gurus will do
    > better than I, though the methodology is very very simple, since all I
    > do is additions, multiplications, and array reads/writes.
    >
    > Have fun with the code.  ;-)
    >


    CORRECTION:

    http://cr.yp.to/primegen.html NOT "primesgen"

    Jabari Zakiya
    jzakiya, Jun 13, 2008
    #2
    1. Advertising

  3. On Jun 13, 1:12 pm, jzakiya <> wrote:

    > The paper presents benchmarks with Ruby 1.9.0-1 (YARV). I would love
    > to see my various prime generators benchmarked with optimized
    > implementations in other languages. I'm hoping Python gurus will do
    > better than I, though the methodology is very very simple, since all I
    > do is additions, multiplications, and array reads/writes.


    After playing a little with it, I managed to get a 32-47% improvement
    on average for the pure Python version, and a 230-650% improvement
    with an extra "import psyco; psyco.full()" (pasted at http://codepad.org/C2nQ8syr)
    The changes are:

    - Replaced range() with xrange()
    - Replaced x**2 with x*x
    - Replaced (a,b) = (c,d) with a=c; b=d
    - Replaced generator expressions with list comprehensions. This was
    the big one for letting psyco do its magic.

    I also tried adding type declarations and running it through Cython
    but the improvement was much less impressive than Psyco. I'm not a
    Pyrex/Cython expert though so I may have missed something obvious.

    George
    George Sakkis, Jun 19, 2008
    #3
  4. jzakiya

    jzakiya Guest

    On Jun 18, 7:58 pm, George Sakkis <> wrote:
    > On Jun 13, 1:12 pm, jzakiya <> wrote:
    >
    > > The paper presents benchmarks with Ruby 1.9.0-1 (YARV).  I would love
    > > to see my variousprimegenerators benchmarked with optimized
    > > implementations in other languages.  I'm hoping Python gurus will do
    > > better than I, though the methodology is very very simple, since all I
    > > do is additions, multiplications, and array reads/writes.

    >
    > After playing a little with it, I managed to get a 32-47% improvement
    > on average for the pure Python version, and a 230-650% improvement
    > with an extra "import psyco; psyco.full()" (pasted athttp://codepad.org/C2nQ8syr)
    > The changes are:
    >
    > - Replaced range() with xrange()
    > - Replaced x**2 with x*x
    > - Replaced (a,b) = (c,d) with a=c; b=d
    > - Replaced generator expressions with list comprehensions. This was
    > the big one for letting psyco do its magic.
    >
    > I also tried adding type declarations and running it through Cython
    > but the improvement was much less impressive than Psyco. I'm not a
    > Pyrex/Cython expert though so I may have missed something obvious.
    >
    > George


    George,

    I took your code and included more efficient/optimized versions of
    SoZ versions P3, P5, P7, and P11.

    I ran the code on my PCLinuxOS, Intel P4, Python 2.4.3 system and
    noted this. The SoZ code run much faster than the SoA in pure Python.
    When psyco is used the SoA is significantly faster than the
    pure Python version. The SoZ versions are faster too, but now
    they are slower than the SoA. You can download the code from


    http://www.4shared.com/dir/7467736/97bd7b71/sharing.html

    It would be interesting to see how this code runs in newer versions
    of Python (Psyco).

    FYI, someone else coded P7 in C on a QuadCore Intel 9650 3.67GHz
    overclocked cpu, using multiple threads, and got it to be faster
    than the SoA, SoE, Here's some of his results (times in seconds).

    Case nPrime7x nPrime7x nPrime7x nPrime7x
    Atkin Zakiya Eratosthenes Zakiya (8 core
    2.5ghz)

    100 billion 52.58 44.27 50.56

    200 b 110.14 92.38 108.99 88.01

    300 b 169.81 140.92 167.47

    400 b 232.34 190.84 228.08 177.72

    500 b 297.44 241.84 291.28

    600 b 364.84 293.92 355.27 273.04

    700 b 434.33 346.97 420.41

    800 b 506.67 400.97 486.72 373.29

    900 b 579.58 456.53 555.09

    1 trillion 654.03 513.11 624.00 479.22


    Jabari
    jzakiya, Jul 19, 2008
    #4
    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:
    4
    Views:
    465
    user923005
    Jun 13, 2008
  2. jzakiya

    Sieve of Zakiya

    jzakiya, Nov 4, 2008, in forum: Python
    Replies:
    7
    Views:
    299
    Mark Dickinson
    Nov 19, 2008
  3. jzakiya
    Replies:
    6
    Views:
    178
    Michael Ulm
    Jun 11, 2008
  4. jzakiya

    Sieve of Zakiya

    jzakiya, Nov 20, 2008, in forum: Ruby
    Replies:
    11
    Views:
    225
  5. Replies:
    0
    Views:
    68
Loading...

Share This Page