Sieve of Zakiya

J

jzakiya

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
 
J

jzakiya

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
 
M

Mark Dickinson


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
 
J

jzakiya

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
 
M

Mark Dickinson

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]])
 
J

jzakiya

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
 
J

jzakiya

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
 
M

Mark Dickinson

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
 

Ask a Question

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

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

Forum statistics

Threads
473,734
Messages
2,569,441
Members
44,832
Latest member
GlennSmall

Latest Threads

Top