Sieve of Zakiya

[jzakiya]

On Nov 18, 6:15 pm, Mark Dickinson <dicki... at gmail.com> wrote:
> On Nov 18, 3:58 pm, jzakiya <jzak... at mail.com> 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[i]]
>
>         # 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[i]:
>             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