"""Eratosthenes.py
Space-efficient version of sieve of Eratosthenes.
D. Eppstein, May 2004.
The main storage of the algorithm is a hash table D with sqrt(n)
nonempty entries for a total of O(sqrt n) space.
At any point in the algorithm, each prime p occupies a cell with key at
most 2n. E.g. by Bertrand's postulate, there is another prime p'
between n/p and 2n/p, and p' can not yet have been included because it
is greater than sqrt n, so key pp' can not be used by any other prime;
therefore p is placed at or before key pp'<2n. Thus, the number of
times p can have been moved from its initial placement at p^2 is < n/p.
The time for the algorithm, up to output n, is O(n) + sum_{prime p <=
sqrt(n)} O(n/p) = O(n log log n). The algorithm also makes a recursive
call, but the recursion only generates primes up to sqrt n so its time
and space is negligible compared to the outer call.
If efficiency is a significant concern it may be better to combine
this idea with segmentation and bitvectors, as in the code by
T. Oliveira e Silva at http://www.ieeta.pt/~tos/software/prime_sieve.html
Thanks to Alex Martelli for the suggestion of keeping one prime
per entry of D, rather than a list of all prime factors of D.
We also include a variant of the sieve that produces a list of all
integers, with their factorizations, and an application of this
variant in the generation of practical numbers.
"""
import unittest
from collections import defaultdict
def primes():
'''Yields the sequence of primes via the Sieve of Eratosthenes.'''
yield 2 # Only even prime. Sieve only odd numbers.
# Generate recursively the sequence of primes up to sqrt(n).
# Each p from the sequence is used to initiate sieving at p*p.
roots = primes()
root = next(roots)
square = root*root
# The main sieving loop.
# We use a hash table D such that D[n]=2p for p a prime factor of n.
# Each prime p up to sqrt(n) appears once as a value in D, and is
# moved to successive odd multiples of p as the sieve progresses.
D = {}
n = 3
while True:
if n >= square: # Time to include another square?
D[square] = root+root
root = next(roots)
square = root*root
if n not in D: # Not witnessed, must be prime.
yield n
else: # Move witness p to next free multiple.
p = D[n]
q = n+p
while q in D:
q += p
del D[n]
D[q] = p
n += 2 # Move on to next odd number.
def FactoredIntegers():
"""
Generate pairs n,F where F is the prime factorization of n.
F is represented as a dictionary in which each prime factor of n
is a key and the exponent of that prime is the corresponding value.
"""
yield 1,{}
i = 2
factorization = defaultdict(dict)
while True:
if i not in factorization: # prime
F = {i:1}
yield i,F
factorization[2*i] = F
elif len(factorization[i]) == 1: # prime power
p,x = next(iter(factorization[i].items()))
F = {p:x+1}
yield i,F
factorization[2*i] = F
factorization[i+p**x][p] = x
del factorization[i]
else:
yield i,factorization[i]
for p,x in factorization[i].items():
q = p**x
iq = i+q
if iq in factorization and p in factorization[iq]:
iq += p**x # skip higher power of p
factorization[iq][p] = x
del factorization[i]
i += 1
def MoebiusSequence():
"""The sequence of values of the Moebius function, OEIS A008683."""
for n,F in FactoredIntegers():
if n > 1 and set(F.values()) != {1}:
yield 0
else:
yield (-1)**len(F)
MoebiusFunctionValues = [None]
MoebiusFunctionIterator = MoebiusSequence()
def MoebiusFunction(n):
"""A functional version of the Moebius sequence.
Efficient only for small values of n."""
while n >= len(MoebiusFunctionValues):
MoebiusFunctionValues.append(next(MoebiusFunctionIterator))
return MoebiusFunctionValues[n]
def isPracticalFactorization(f):
"""Test whether f is the factorization of a practical number."""
f = sorted(f.items())
sigma = 1
for p,x in f:
if sigma < p - 1:
return False
sigma *= (p**(x+1)-1)//(p-1)
return True
def PracticalNumbers():
"""Generate the sequence of practical (or panarithmic) numbers."""
for x,f in FactoredIntegers():
if isPracticalFactorization(f):
yield x
# If run standalone, perform unit tests
class SieveTest(unittest.TestCase):
def testPrime(self):
"""Test that the first few primes are generated correctly."""
G = primes()
for p in [2,3,5,7,11,13,17,19,23,29,31,37]:
self.assertEqual(p,next(G))
def testPractical(self):
"""Test that the first few practical nos are generated correctly."""
G = PracticalNumbers()
for p in [1,2,4,6,8,12,16,18,20,24,28,30,32,36]:
self.assertEqual(p,next(G))
if __name__ == "__main__":
unittest.main()