This algorithm enumerates the twin-primes without performing explicit primality tests.
A twin-prime witness is an integer w such that both 6w−1 and 6w+1 are both prime.
A characterization discussed in OEIS A002822 and in work of Francesca Balestrieri is that w is a twin-prime witness if and only if it cannot be expressed in any of the forms
The algorithm systematically enumerates all values of the form 6ab±a±b. It maintains a priority queue of such values and emits the integers that are not covered. These uncovered integers are precisely the twin-prime witnesses, from which the corresponding twin-prime pairs 6w−1,6w+1 are produced.
import heapq
class TwinPrimeWalk:
def __init__(self):
pass
def encode_sigma(self, c, d, sigma_c, sigma_d):
n = c * d
f = 6*n+sigma_c*c+sigma_d*d
t = (sigma_c+1)//2+sigma_d+1
return (f, -n, c, d, t)
def encode(self, c,d, t):
return self.encode_sigma(c,d, (t%2)*2-1, (t//2)*2-1)
def twin_primes(self):
yield 3
q = []
w = 0
last_sqn=1
heapq.heappush(q, self.encode(1,1,0))
while len(q) > 0:
r = heapq.heappop(q)
f, _n, c, d, t = r
n=-_n
for ww in range(w+1, f):
yield(6*ww-1)
yield(6*ww+1)
w = f
if t == 0:
heapq.heappush(q, self.encode(c, d, 1))
heapq.heappush(q, self.encode(c, d, 2))
heapq.heappush(q, self.encode(c, d, 3))
heapq.heappush(q, self.encode(c, d+1, 0))
while n > last_sqn**2:
sqn = last_sqn+1
heapq.heappush(q, self.encode(sqn, sqn, 0))
last_sqn = sqn
[ w for i, w in zip(range(0, 100), TwinPrimeWalk().twin_primes()) ]
[1] OEIS A002822 - the OEIS sequence that is the set of witnesses of twin primes
[2] F. Balestrieri, An Equivalent Problem To The Twin Prime Conjecture, arXiv:1106.6050v1 [math.GM], 2011.
[3] J. Seymour, "The Sieve of Balestrieri", a visualisation
[4] Suzuki, M. (2000). Alternative formulations of the twin prime problem. The American Mathematical Monthly, 107(1), 55-56. (h/t u/davidjohnpaul for finding this)