On Nov 18, 3:58 pm, jzakiya <[EMAIL PROTECTED]> 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]])
--
http://mail.python.org/mailman/listinfo/python-list
