from sympy import primerange from heapq import heappush, heappop top = 10**9 heap = [(3, 3), (4, 2)] last = 2 primes = primerange(5, top) seen = {} while heap: n, p = heappop(heap) gap = n - last if not (gap in seen): seen[gap] = last last = n m = n * p if m < top: heappush(heap, (m, p)) if n == p: try: m = primes.next() if m < top: heappush(heap, (m, m)) except: pass for x, y in sorted(seen.items()): print "a(%d) = %d" % (x, y)