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:
] = [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[:
] = [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]])
