Sieve of Zakiya

Discussion in 'Python' started by jzakiya, Nov 4, 2008.

  1. jzakiya

    jzakiya Guest

    Update: 2008/11/03

    Architecture & coding improvements. Renamed generators.

    I am 90% finished writing up a mathematical analysis of my method.
    In the process I found an architectural optimization to the sieve
    process which is incorporated in the new code. Complexity analysis
    showing other interesting stuff for each generator.

    When I finish I will post paper here with the code:

    www.4shared.com/account/dir/7467736/97bd7b71/sharing

    Jabari
     
    jzakiya, Nov 4, 2008
    #1
    1. Advertising

  2. jzakiya

    jzakiya Guest

    On Nov 4, 4:12 pm, jzakiya <> wrote:
    > Update: 2008/11/03
    >
    > Architecture & coding improvements. Renamed generators.
    >
    > I am 90% finished writing up a mathematical analysis of my method.
    > In the process I found an architectural optimization to the sieve
    > process which is incorporated in the new code. Complexity analysis
    > showing other interesting stuff for each generator.
    >
    > When I finish I will post paper here with the code:
    >
    > www.4shared.com/account/dir/7467736/97bd7b71/sharing
    >
    > Jabari


    Another update to SoZ code.

    Fixed coding "oversight" which caused the SoZ code running under psyco
    to be slower than SoA, though SoZ was faster than SoA without psyco.
    Now SoZ is consistently faster by an appreciable amount either way.

    Some other coding enhancements to optimize algorithmic performance.

    Enjoy! :)

    Jabari
     
    jzakiya, Nov 18, 2008
    #2
    1. Advertising

  3. On Nov 18, 7:53 am, jzakiya <> wrote:
    > www.4shared.com/account/dir/7467736/97bd7b71/sharing


    From the introduction to the paper:

    "Thus began a process that culminated in my developing a new class
    of Number Theory Sieves (NTS) to generate prime numbers, and test
    primality of numbers, that use minimum memory, are simple to code,
    and are much faster than all previously known methods."

    That's quite a claim! You might consider toning this down
    a little if you want to be taken seriously. :)

    How does your prime-generating algorithm stack up against
    Bernstein's primegen package?

    http://cr.yp.to/primegen.html

    Mark
     
    Mark Dickinson, Nov 18, 2008
    #3
  4. jzakiya

    jzakiya Guest

    On Nov 18, 5:01 am, Mark Dickinson <> wrote:
    > On Nov 18, 7:53 am, jzakiya <> wrote:
    >
    > >www.4shared.com/account/dir/7467736/97bd7b71/sharing

    >
    > From the introduction to the paper:
    >
    > "Thus began a process that culminated in my developing a new class
    > of Number Theory Sieves (NTS) to generate prime numbers, and test
    > primality of numbers, that use minimum memory, are simple to code,
    > and are much faster than all previously known methods."
    >
    > That's quite a claim!  You might consider toning this down
    > a little if you want to be taken seriously.  :)
    >
    > How does your prime-generating algorithm stack up against
    > Bernstein's primegen package?
    >
    > http://cr.yp.to/primegen.html
    >
    > Mark


    Hi Mark,

    Someone has done a multi-core C implementation (for AMD-X2 and Intel
    quad and 8-core systems using Intel and GCC compilers) and the SoZ
    prime generators are faster than the Sieve of Atkin (SoA) and Sieve of
    Eratosthenes (SoE).

    I think(?) he tried to run Bernstein's primegen, but its old code
    written before multi-core and widely used 64-bit systems, so he wrote
    his own multi-core/threaded version for his systems. The SoZs again,
    are faster.

    I've implemented the SoZ generators in Forth, Ruby, and Python, and
    all show the same results, to be faster than the SoA and SoE in those
    languages, along with C.

    Now that I've corrected the Python code (I am NOT a native Python
    programmer) it now beats the SoA implementation included in the code
    using psyco too (for 2.4.3, I don't have it for other versions).

    Run the code and see for yourself! :)

    I am writing another paper explaining some of the mathematical basis
    for the SoZ, with complexity analysis, but I keep finding
    "interesting" features about the underlying math, which I hope real
    mathematicians will investigate and reveal what's going on here.
    I want to release at least a first version before December.

    But the proof is in the pudding, and the results speak for themselves.
    It's an amazingly simple algorithm, so tear it apart.

    Jabari
     
    jzakiya, Nov 18, 2008
    #4
  5. On Nov 18, 3:58 pm, jzakiya <> wrote:
    > I am writing another paper explaining some of the mathematical basis
    > for the SoZ, with complexity analysis, but I keep finding
    > "interesting" features about the underlying math, which I hope real
    > mathematicians will investigate and reveal what's going on here.


    Well, what's going on here is that you've rediscovered the idea
    of a wheel. I don't know who first came up with the idea of
    speeding up the sieve of Eratosthenes with a wheel, but it's
    certainly been around for a while. So unfortunately your
    ideas, while interesting, aren't new. Don't be put off,
    though: finding that your amazing new discovery is already
    well known is a pretty common activity when doing mathematics. :)

    A good reference for this and related topics is the book
    "Prime Numbers: A Computational Perspective", by Crandall
    and Pomerance. See Chapter 3, especially section 3.1.

    > Run the code and see for yourself! :)


    Using Python to measure algorithm performance is kinda
    crazy, but since you insist, here's a version of the sieve
    of Eratosthenes that I wrote a while back. On my machine,
    the function wheelSieve comes out around 60% faster than
    your 'SoZP11', for inputs of around 10**6 or so.

    Mark


    from math import sqrt
    from bisect import bisect_left

    def basicSieve(n):
    """Given a positive integer n, generate the primes < n."""
    s = [1]*n
    for p in xrange(2, 1+int(sqrt(n-1))):
    if s[p]:
    a = p*p
    s[a::p] = [0] * -((a-n)//p)
    for p in xrange(2, n):
    if s[p]:
    yield p

    class Wheel(object):
    """Class representing a wheel.

    Attributes:
    primelimit -> wheel covers primes < primelimit.
    For example, given a primelimit of 6
    the wheel primes are 2, 3, and 5.
    primes -> list of primes less than primelimit
    size -> product of the primes in primes; the modulus of the
    wheel
    units -> list of units modulo size
    phi -> number of units

    """
    def __init__(self, primelimit):
    self.primelimit = primelimit
    self.primes = list(basicSieve(primelimit))

    # compute the size of the wheel
    size = 1
    for p in self.primes:
    size *= p
    self.size = size

    # find the units by sieving
    units = [1] * self.size
    for p in self.primes:
    units[::p] = [0]*(self.size//p)
    self.units = [i for i in xrange(self.size) if units]

    # number of units
    self.phi = len(self.units)

    def to_index(self, n):
    """Compute alpha(n), where alpha is an order preserving map
    from the set of units modulo self.size to the nonnegative
    integers.

    If n is not a unit, the index of the first unit greater than n
    is given."""
    return bisect_left(self.units, n % self.size) + n // self.size
    * self.phi

    def from_index(self, i):
    """Inverse of to_index."""

    return self.units[i % self.phi] + i // self.phi * self.size

    def wheelSieveInner(n, wheel):
    """Given a positive integer n and a wheel, find the wheel indices
    of
    all primes that are less than n, and that are units modulo the
    wheel modulus.
    """

    # renaming to avoid the overhead of attribute lookups
    U = wheel.units
    wS = wheel.size
    # inverse of unit map
    UI = dict((u, i) for i, u in enumerate(U))
    nU = len(wheel.units)

    sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n

    # corresponding index (index of next unit up)
    sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
    upper = bisect_left(U, n % wS) + n//wS*nU
    ind2 = bisect_left(U, 2 % wS) + 2//wS*nU

    s = [1]*upper
    for i in xrange(ind2, sqrti):
    if s:
    q = i//nU
    u = U[i%nU]
    p = q*wS+u
    u2 = u*u
    aq, au = (p+u)*q+u2//wS, u2%wS
    wp = p * nU
    for v in U:
    # eliminate entries corresponding to integers
    # congruent to p*v modulo p*wS
    uvr = u*v%wS
    m = aq + (au > uvr)
    bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
    s[bot::wp] = [0]*-((bot-upper)//wp)
    return s

    def wheelSieve(n, wheel=Wheel(10)):
    """Given a positive integer n, generate the list of primes <=
    n."""
    n += 1
    wS = wheel.size
    U = wheel.units
    nU = len(wheel.units)
    s = wheelSieveInner(n, wheel)
    return (wheel.primes[:bisect_left(wheel.primes, n)] +
    [p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
    + 2//wS*nU, len(s)) if s[p]])
     
    Mark Dickinson, Nov 18, 2008
    #5
  6. jzakiya

    jzakiya Guest

    On Nov 18, 6:15 pm, Mark Dickinson <> wrote:
    > On Nov 18, 3:58 pm, jzakiya <> wrote:
    >
    > > I am writing another paper explaining some of the mathematical basis
    > > for the SoZ, with complexity analysis, but I keep finding
    > > "interesting" features about the underlying math, which I hope real
    > > mathematicians will investigate and reveal what's going on here.

    >
    > Well, what's going on here is that you've rediscovered the idea
    > of a wheel.  I don't know who first came up with the idea of
    > speeding up the sieve of Eratosthenes with a wheel, but it's
    > certainly been around for a while.  So unfortunately your
    > ideas, while interesting, aren't new.  Don't be put off,
    > though: finding that your amazing new discovery is already
    > well known is a pretty common activity when doing mathematics. :)
    >
    > A good reference for this and related topics is the book
    > "Prime Numbers: A Computational Perspective", by Crandall
    > and Pomerance.  See Chapter 3, especially section 3.1.
    >
    > > Run the code and see for yourself! :)

    >
    > Using Python to measure algorithm performance is kinda
    > crazy, but since you insist, here's a version of the sieve
    > of Eratosthenes that I wrote a while back.  On my machine,
    > the function wheelSieve comes out around 60% faster than
    > your 'SoZP11', for inputs of around 10**6 or so.
    >
    > Mark
    >
    > from math import sqrt
    > from bisect import bisect_left
    >
    > def basicSieve(n):
    >     """Given a positive integer n, generate the primes < n."""
    >     s = [1]*n
    >     for p in xrange(2, 1+int(sqrt(n-1))):
    >         if s[p]:
    >             a = p*p
    >             s[a::p] = [0] * -((a-n)//p)
    >     for p in xrange(2, n):
    >         if s[p]:
    >             yield p
    >
    > class Wheel(object):
    >     """Class representing a wheel.
    >
    >     Attributes:
    >        primelimit -> wheel covers primes < primelimit.
    >        For example, given a primelimit of 6
    >        the wheel primes are 2, 3, and 5.
    >        primes -> list of primes less than primelimit
    >        size -> product of the primes in primes;  the modulus of the
    > wheel
    >        units -> list of units modulo size
    >        phi -> number of units
    >
    >     """
    >     def __init__(self, primelimit):
    >         self.primelimit = primelimit
    >         self.primes = list(basicSieve(primelimit))
    >
    >         # compute the size of the wheel
    >         size = 1
    >         for p in self.primes:
    >             size *= p
    >         self.size = size
    >
    >         # find the units by sieving
    >         units = [1] * self.size
    >         for p in self.primes:
    >             units[::p] = [0]*(self.size//p)
    >         self.units = [i for i in xrange(self.size) if units]
    >
    >         # number of units
    >         self.phi = len(self.units)
    >
    >     def to_index(self, n):
    >         """Compute alpha(n), where alpha is an order preserving map
    >         from the set of units modulo self.size to the nonnegative
    > integers.
    >
    >         If n is not a unit, the index of the first unit greater than n
    >         is given."""
    >         return bisect_left(self.units, n % self.size) + n // self..size
    > * self.phi
    >
    >     def from_index(self, i):
    >         """Inverse of to_index."""
    >
    >         return self.units[i % self.phi] + i // self.phi * self.size
    >
    > def wheelSieveInner(n, wheel):
    >     """Given a positive integer n and a wheel, find the wheel indices
    > of
    >     all primes that are less than n, and that are units modulo the
    >     wheel modulus.
    >     """
    >
    >     # renaming to avoid the overhead of attribute lookups
    >     U = wheel.units
    >     wS = wheel.size
    >     # inverse of unit map
    >     UI = dict((u, i) for i, u in enumerate(U))
    >     nU = len(wheel.units)
    >
    >     sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n
    >
    >     # corresponding index (index of next unit up)
    >     sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
    >     upper = bisect_left(U, n % wS) + n//wS*nU
    >     ind2 = bisect_left(U, 2 % wS) + 2//wS*nU
    >
    >     s = [1]*upper
    >     for i in xrange(ind2, sqrti):
    >         if s:
    >             q = i//nU
    >             u = U[i%nU]
    >             p = q*wS+u
    >             u2 = u*u
    >             aq, au = (p+u)*q+u2//wS, u2%wS
    >             wp = p * nU
    >             for v in U:
    >                 # eliminate entries corresponding to integers
    >                 # congruent to p*v modulo p*wS
    >                 uvr = u*v%wS
    >                 m = aq + (au > uvr)
    >                 bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
    >                 s[bot::wp] = [0]*-((bot-upper)//wp)
    >     return s
    >
    > def wheelSieve(n, wheel=Wheel(10)):
    >     """Given a positive integer n, generate the list of primes <=
    > n."""
    >     n += 1
    >     wS = wheel.size
    >     U = wheel.units
    >     nU = len(wheel.units)
    >     s = wheelSieveInner(n, wheel)
    >     return (wheel.primes[:bisect_left(wheel.primes, n)] +
    >             [p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
    >              + 2//wS*nU, len(s)) if s[p]])


    This is not an apples-to-apples comparison. You are using an SoE to
    create a composite wheel sieve. Well, you could do the same with a SoZ
    generator too.

    According to my research, the Sieve of Atkin (SoA) was recognized to
    be the fastest contemporary prime sieve method. My SoZ clearly beats
    it. By inference then, you are saying your method is also faster than
    the SoA. Is that what you are saying?

    I googled (again) 'wheel sieves' and this gave a good explanation:

    http://primes.utm.edu/glossary/page.php?sort=WheelFactorization

    Some of what it explains has some of my properties, but I have not
    seen any method that eliminates non-primes as I do, and certainly not
    faster than the SoA.

    It would be very helpful if you inserted your wheel sieve into my
    Sieves.py file and created a comparison that can be run like I
    provided. You didn't present any real numbers, and a reproducible
    benchmark to see how you got your results. Also, 1M (10^6) is kind of
    small. My SoZs start to really shine when you get into the 1B (10^9)
    range, and higher. What are your results in this range?

    My SoZ mathematical analysis shows why its algorithmically faster than
    both the SoA & SoE, to explain its demonstrated computational
    superiority. Again, hopefully professional mathematicians will become
    interested in some of the features and characteristics I've found in
    the method's elements, and do a better and more thorough analysis of
    what's going on than I'm able to do presently.

    Jabari
     
    jzakiya, Nov 19, 2008
    #6
  7. jzakiya

    jzakiya Guest

    On Nov 18, 6:15 pm, Mark Dickinson <> wrote:
    > On Nov 18, 3:58 pm, jzakiya <> wrote:
    >
    > > I am writing another paper explaining some of the mathematical basis
    > > for the SoZ, with complexity analysis, but I keep finding
    > > "interesting" features about the underlying math, which I hope real
    > > mathematicians will investigate and reveal what's going on here.

    >
    > Well, what's going on here is that you've rediscovered the idea
    > of a wheel.  I don't know who first came up with the idea of
    > speeding up the sieve of Eratosthenes with a wheel, but it's
    > certainly been around for a while.  So unfortunately your
    > ideas, while interesting, aren't new.  Don't be put off,
    > though: finding that your amazing new discovery is already
    > well known is a pretty common activity when doing mathematics. :)
    >
    > A good reference for this and related topics is the book
    > "Prime Numbers: A Computational Perspective", by Crandall
    > and Pomerance.  See Chapter 3, especially section 3.1.
    >
    > > Run the code and see for yourself! :)

    >
    > Using Python to measure algorithm performance is kinda
    > crazy, but since you insist, here's a version of the sieve
    > of Eratosthenes that I wrote a while back.  On my machine,
    > the function wheelSieve comes out around 60% faster than
    > your 'SoZP11', for inputs of around 10**6 or so.
    >
    > Mark
    >
    > from math import sqrt
    > from bisect import bisect_left
    >
    > def basicSieve(n):
    >     """Given a positive integer n, generate the primes < n."""
    >     s = [1]*n
    >     for p in xrange(2, 1+int(sqrt(n-1))):
    >         if s[p]:
    >             a = p*p
    >             s[a::p] = [0] * -((a-n)//p)
    >     for p in xrange(2, n):
    >         if s[p]:
    >             yield p
    >
    > class Wheel(object):
    >     """Class representing a wheel.
    >
    >     Attributes:
    >        primelimit -> wheel covers primes < primelimit.
    >        For example, given a primelimit of 6
    >        the wheel primes are 2, 3, and 5.
    >        primes -> list of primes less than primelimit
    >        size -> product of the primes in primes;  the modulus of the
    > wheel
    >        units -> list of units modulo size
    >        phi -> number of units
    >
    >     """
    >     def __init__(self, primelimit):
    >         self.primelimit = primelimit
    >         self.primes = list(basicSieve(primelimit))
    >
    >         # compute the size of the wheel
    >         size = 1
    >         for p in self.primes:
    >             size *= p
    >         self.size = size
    >
    >         # find the units by sieving
    >         units = [1] * self.size
    >         for p in self.primes:
    >             units[::p] = [0]*(self.size//p)
    >         self.units = [i for i in xrange(self.size) if units]
    >
    >         # number of units
    >         self.phi = len(self.units)
    >
    >     def to_index(self, n):
    >         """Compute alpha(n), where alpha is an order preserving map
    >         from the set of units modulo self.size to the nonnegative
    > integers.
    >
    >         If n is not a unit, the index of the first unit greater than n
    >         is given."""
    >         return bisect_left(self.units, n % self.size) + n // self..size
    > * self.phi
    >
    >     def from_index(self, i):
    >         """Inverse of to_index."""
    >
    >         return self.units[i % self.phi] + i // self.phi * self.size
    >
    > def wheelSieveInner(n, wheel):
    >     """Given a positive integer n and a wheel, find the wheel indices
    > of
    >     all primes that are less than n, and that are units modulo the
    >     wheel modulus.
    >     """
    >
    >     # renaming to avoid the overhead of attribute lookups
    >     U = wheel.units
    >     wS = wheel.size
    >     # inverse of unit map
    >     UI = dict((u, i) for i, u in enumerate(U))
    >     nU = len(wheel.units)
    >
    >     sqroot = 1+int(sqrt(n-1)) # ceiling of square root of n
    >
    >     # corresponding index (index of next unit up)
    >     sqrti = bisect_left(U, sqroot % wS) + sqroot//wS*nU
    >     upper = bisect_left(U, n % wS) + n//wS*nU
    >     ind2 = bisect_left(U, 2 % wS) + 2//wS*nU
    >
    >     s = [1]*upper
    >     for i in xrange(ind2, sqrti):
    >         if s:
    >             q = i//nU
    >             u = U[i%nU]
    >             p = q*wS+u
    >             u2 = u*u
    >             aq, au = (p+u)*q+u2//wS, u2%wS
    >             wp = p * nU
    >             for v in U:
    >                 # eliminate entries corresponding to integers
    >                 # congruent to p*v modulo p*wS
    >                 uvr = u*v%wS
    >                 m = aq + (au > uvr)
    >                 bot = (m + (q*v + u*v//wS - m) % p) * nU + UI[uvr]
    >                 s[bot::wp] = [0]*-((bot-upper)//wp)
    >     return s
    >
    > def wheelSieve(n, wheel=Wheel(10)):
    >     """Given a positive integer n, generate the list of primes <=
    > n."""
    >     n += 1
    >     wS = wheel.size
    >     U = wheel.units
    >     nU = len(wheel.units)
    >     s = wheelSieveInner(n, wheel)
    >     return (wheel.primes[:bisect_left(wheel.primes, n)] +
    >             [p//nU*wS + U[p%nU] for p in xrange(bisect_left(U, 2 % wS)
    >              + 2//wS*nU, len(s)) if s[p]])


    This is not an apples-to-apples comparison. You are using an SoE to
    create a composite wheel sieve. Well, you could do the same with a SoZ
    generator too.

    According to my research, the Sieve of Atkin (SoA) was recognized to
    be the fastest contemporary prime sieve method. My SoZ clearly beats
    it. By inference then, you are saying your method is also faster than
    the SoA. Is that what you are saying?

    I googled (again) 'wheel sieves' and this gave a good explanation:

    http://primes.utm.edu/glossary/page.php?sort=WheelFactorization

    Some of what it explains has some of my properties, but I have not
    seen any method that eliminates non-primes as I do, and certainly not
    faster than the SoA.

    It would be very helpful if you inserted your wheel sieve into my
    Sieves.py file and created a comparison that can be run like I
    provided. You didn't present any real numbers, and a reproducible
    benchmark to see how you got your results. Also, 1M (10^6) is kind of
    small. My SoZs start to really shine when you get into the 1B (10^9)
    range, and higher. What are your results in this range?

    My SoZ mathematical analysis shows why its algorithmically faster than
    both the SoA & SoE, to explain its demonstrated computational
    superiority. Again, hopefully professional mathematicians will become
    interested in some of the features and characteristics I've found in
    the method's elements, and do a better and more thorough analysis of
    what's going on than I'm able to do presently.

    Jabari
     
    jzakiya, Nov 19, 2008
    #7
  8. On Nov 19, 5:16 am, jzakiya <> wrote:
    [Lots of stuff snipped]
    > My SoZ mathematical analysis shows why its algorithmically faster than
    > both the SoA & SoE, to explain its demonstrated computational
    > superiority. Again, hopefully professional mathematicians will become
    > interested in some of the features and characteristics I've found in
    > the method's elements, and do a better and more thorough analysis of
    > what's going on than I'm able to do presently.


    This is getting off-topic for comp.lang.python. I suggest you
    take it to sci.math.

    Mark
     
    Mark Dickinson, Nov 19, 2008
    #8
    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:
    577
    jzakiya
    Jul 19, 2008
  2. jzakiya
    Replies:
    4
    Views:
    500
    user923005
    Jun 13, 2008
  3. jzakiya
    Replies:
    6
    Views:
    209
    Michael Ulm
    Jun 11, 2008
  4. jzakiya

    Sieve of Zakiya

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

Share This Page