Personal Web Page

Thursday, July 30, 2009

 

Lazy Primes Sieve in Python

I read this interesting paper, The Genuine Sieve of Eratosthenes (PDF) by Melissa O'Neill, when it was discussed on Lambda the Ultimate. I came back to look at it again today because I wanted to improve my Haskell fluency, which I only seem to exercise when I read papers like this. I decided to try to correctly construe the Haskell code in the paper by implementing the algorithms in Python. I focused on reproducing the laziness but not the functional style because I could substitute iterators for lazy lists but had to use conventional control structures for all the pattern matching and tail calls.


#!/usr/bin/env python
# Author: Jared Brothers
#
# A Python version of the Haskell code from
# "The Genuine Sieve of Eratosthenes"
# www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
#
import heapq, itertools, operator

def sieve(xs):
"""Generate the prime numbers, given an iterable of candidate numbers.
Cross off multiples of prime numbers incrementally using iterators."""
table = []
while 1:
x = xs.next()
if table == [] or x < table[0][0]:
yield x
xs, ys = itertools.tee(xs)
timesx = (lambda x: lambda y: x*y)(x)
heapq.heappush(table, (x**2, itertools.imap(timesx, ys)))
else:
while table[0][0] <= x:
heapq.heapreplace(table, (table[0][1].next(), table[0][1]))

def wheel(factors=[2, 3, 5, 7], next=11):
"""Generate the distances between numbers not divisible by a list of small
primes, from the next prime up to the product of the list."""
circumference = reduce(operator.mul, factors)
prev = next
next += 1
end = next + circumference
while next < end:
if not any(next % factor == 0 for factor in factors):
yield next - prev
prev = next
next += 1

def spin(factors=[2, 3, 5, 7], next=11):
"""Generate candidates by making a wheel and cycling through it."""
for gap in itertools.cycle(wheel(factors, next)):
yield next
next += gap

def primes(k=5):
"""Generate primes with the sieve and wheel factorization, which filters
multiples of the first k primes."""
smallprimes = list(itertools.islice(sieve(itertools.count(2)), k + 1))
factors = smallprimes[:-1]
next = smallprimes[-1]
return itertools.chain(factors, sieve(spin(factors, next)))

Labels: , , ,


Sunday, January 4, 2009

 

Domino counting problem

My brother is writing a combinatorics paper about playing with dominoes. Pure mathematics is a lot of fun that way, but with more Hamiltonian paths of graphs embedded on a torus. Here is a Python script I wrote for him last night to solve a counting problem.


#!/usr/bin/env python
# Author: Jared Brothers
# Find the orderings of a set of dominoes such that
# adjacent dominoes share a number and the other two
# numbers differ by one.

def search(start):
""" Depth first search, generating all solutions. """
fringe = [start]
while fringe:
s = fringe.pop()
if s.check():
yield s
fringe.extend(s.successors())

class state():
""" A state in the search space. """
def __init__(self, d, r):
""" The ordered dominoes are in done (list),
and the unordered dominoes are in rest (set). """
self.done = d
self.rest = r

def __str__(self):
return str(self.done)

def check(self):
""" Is this a solution? """
return not self.rest

def successors(self):
""" The states you can get to by moving a domino from rest to done. """
if self.done:
# Try only dominoes adjacent to the previous one,
a,b = self.done[-1]
a1 = (a + 1) % n
a2 = (a - 1) % n
b1 = (b + 1) % n
b2 = (b - 1) % n
adjs = set([(a,b1), (b1,a), (a,b2), (b2,a), (b,a1), (a1,b), (b,a2), (a2,b)])
for x in self.rest.intersection(adjs):
d = self.done[:]
d.append(x)
r = self.rest.copy()
r.discard(x)
yield state(d, r)
else:
# or all the rest if nothing has been done yet.
for x in sorted(self.rest):
d = [x]
r = self.rest.copy()
r.discard(x)
yield state(d, r)

if __name__ == '__main__':
n = 7 # The number of numbers. The number of dominoes is sum(range(n + 1)).
start = state([], set((a, b) for a in range(n) for b in range(a + 1)))
for s in search(start):
print(s)

Labels: ,


My Photo
Name: Jared Brothers
Location: New York, NY, United States

Links

Archives

November 2008   December 2008   January 2009   June 2009   July 2009  

This page is powered by Blogger. Isn't yours?

Subscribe to Posts [Atom]