## Prime Chains

### July 4, 2017

Our solution uses breadth-first search through the neighborhood of the starting prime, continuing until the target prime is found. First we write a simple function for checking primality:

```(define (prime? n)
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((f 2) (w 0))
(if (< n (* f f)) #t
(if (zero? (modulo n f)) #f
(loop (+ f (vector-ref wheel w))
(if (= w 10) 3 (+ w 1))))))))```

That performs trial division using a 2,3,5-wheel; it is sufficient for five-digit numbers, as in the sample, but you’ll want something faster if you use larger numbers. The neighbors of a number are those that can be reached by altering a single digit; a k-digit number has 9k − 1 neighbors:

```(define (neighbors n)
(let ((ds (digits n)))
(let outer ((prefix (list)) (digit (car ds)) (suffix (cdr ds)) (ns (list)))
(if (null? suffix)
(let inner ((rs (remove digit (range (if (pair? prefix) 0 1) 10))) (xs (list)))
(if (null? rs) (append xs ns)
(inner (cdr rs) (cons (undigits (append prefix (list (car rs)))) xs))))
(let inner ((rs (remove digit (range (if (pair? prefix) 0 1) 10))) (xs (list)))
(if (null? rs) (outer (append prefix (list digit)) (car suffix) (cdr suffix) (append xs ns))
(inner (cdr rs) (cons (undigits (append prefix (list (car rs)) suffix)) xs))))))))```

For example, the neighbors of 71549, and the prime neighbors of 71549, can be computed as shown below:

```> (neighbors 71549)
(71548 71547 71546 71545 71544 71543 71542 71541 71540 71599
71589 71579 71569 71559 71539 71529 71519 71509 71949 71849
71749 71649 71449 71349 71249 71149 71049 79549 78549 77549
76549 75549 74549 73549 72549 70549 91549 81549 61549 51549
41549 31549 21549 11549)
> (filter prime? (neighbors 71549))
(71569 71849 71249 79549 77549 70549 41549 11549)```

At each step in the breadth-first search, we take a list of chains (a list of lists, each containing the starting prime at its tail end) and extend each chain in the list by its prime neighbors:

```(define (extend chains)
(let loop ((chains chains) (xss (list)))
(if (null? chains) xss
(loop (cdr chains)
(append (map (lambda (n) (cons n (car chains)))
(filter prime? (neighbors (caar chains))))
xss)))))```

The list of chains quickly grows large:

```> (extend '((71549)))
((71569 71549) (71849 71549) (71249 71549) (79549 71549)
(77549 71549) (70549 71549) (41549 71549) (11549 71549))
> (extend (extend '((71549))))
((11579 11549 71549) (11519 11549 71549) (11149 11549 71549)
(14549 11549 71549) (71549 11549 71549) (41549 11549 71549)
(41543 41549 71549) (41579 41549 71549) (41539 41549 71549)
(41519 41549 71549) (41849 41549 71549) (41149 41549 71549)
(49549 41549 71549) (46549 41549 71549) (44549 41549 71549)
(71549 41549 71549) (11549 41549 71549) (70589 70549 71549)
(70529 70549 71549) (70949 70549 71549) (70849 70549 71549)
(70249 70549 71549) (79549 70549 71549) (77549 70549 71549)
(71549 70549 71549) (50549 70549 71549) (20549 70549 71549)
(77543 77549 71549) (77569 77549 71549) (77509 77549 71549)
(77849 77549 71549) (77249 77549 71549) (79549 77549 71549)
(71549 77549 71549) (70549 77549 71549) (97549 77549 71549)
(37549 77549 71549) (79589 79549 71549) (79579 79549 71549)
(79559 79549 71549) (79349 79549 71549) (77549 79549 71549)
(71549 79549 71549) (70549 79549 71549) (49549 79549 71549)
(71209 71249 71549) (71849 71249 71549) (71549 71249 71549)
(77249 71249 71549) (76249 71249 71549) (70249 71249 71549)
(91249 71249 71549) (31249 71249 71549) (71843 71849 71549)
(71899 71849 71549) (71879 71849 71549) (71809 71849 71549)
(71549 71849 71549) (71249 71849 71549) (77849 71849 71549)
(73849 71849 71549) (70849 71849 71549) (41849 71849 71549)
(31849 71849 71549) (71563 71569 71549) (71549 71569 71549)
(71069 71569 71549) (78569 71569 71549) (77569 71569 71549)
(81569 71569 71549) (21569 71569 71549))```

Finally, we recursively extend the chains until we find one or more chains headed by the target prime:

```(define (chains start target)
(let loop ((chains (list (list start))))
(let ((xss (extend chains)))
(let ((zss (filter (lambda (xs) (= (car xs) target)) xss)))
(if (pair? zss) (map reverse zss) (loop xss))))))```

Here are two examples:

```> (chains 71549 10067)
((71549 71569 71069 11069 10069 10067))
> (chains 71549 18773)
((71549 11549 11519 11719 18719 18713 18773)
(71549 11549 11149 18149 18143 18743 18773)
(71549 11549 11149 18149 18749 18743 18773))```

We use range, digits and undigits from the Standard Prelude. You can run the program at http://ideone.com/0nirN1.

Pages: 1 2

### 2 Responses to “Prime Chains”

1. Paul said

In Python. This is taken from the book Python Algorithms by Hetland. In the book code is given for finding chains of words from a dictionary. I only had to change the method “variants”. The A* method is used to find the shortest path with a distance heuristic that counts the number of mismatches in the numbers.

```def a_star(G, s, t, h):
P, Q = {}, [(h(s), None, s)]                # Pred and queue w/heuristic
while Q:                                    # Still unprocessed nodes?
d, p, u = heappop(Q)                    # Node with lowest heuristic
if u in P: continue                     # Already visited? Skip it
P[u] = p                                # Set path predecessor
if u == t: return d - h(t), P           # Arrived! Ret. dist and preds
for v in G[u]:                          # Go through all neighbors
w = G[u][v] - h(u) + h(v)           # Modify weight wrt heuristic
heappush(Q, (d + w, u, v))          # Add to queue, w/heur as pri
return inf, None                            # Didn't get to t

class PrimeChain(object):

def __init__(self):
self.M = {}

def variants(self, strn):
va = list(strn)
for i, c in enumerate(va):
digits = "123456789" if i == 0 else "0123456789"
for d in digits:
if d == c:
continue
va[i] = d
newstrn = "".join(va)
if is_prime(int(newstrn)):
yield newstrn
va[i] = c

def __getitem__(self, strn):
if strn not in self.M:
self.M[strn] = dict.fromkeys(self.variants(strn), 1)
return self.M[strn]

def heuristic(self, u, v):
return sum(a != b for a, b in zip(u, v))

def chain(self, s, t, h=None):
if h is None:
def h(v):
return self.heuristic(v, t)
_, P = a_star(self, s, t, h)
if P is None:
return [s, None, t]
u, p = t, []
while u is not None:
p.append(u)
u = P[u]
p. reverse()
return p

G = PrimeChain()
print(G.chain("71549", "10067"))
p1 = 76520965667
p2 = 59340926171
print(G.chain(str(p1), str(p2)))

"""
['71549', '71569', '71069', '11069', '10069', '10067']
['76520965667', '70520965667', '70540965667', '79540965667', '79040965667',
'79040965657', '79040965651', '79040966651', '79040966671', '59040966671',
'59040926671', '590409261 71', '59340926171']
"""
```
2. Rutger said

Hi,

Created a Depth-First search algorithm. This, off-course, goes bananas without a proper heuristic on getting closer to the target prime. This was done using the Hamming distance. Also set a limit of max depth (best) = 100 due to recursion depth limits.

Also found out that the is_prime function has many calls with the same argument, so cached the results using memoization.

I recon BFS would be a better approach judging from your solutions, anyway here’s my result:

```import random
import string

def memoize(f):
memo = {}
def helper(x):
if x not in memo:
memo[x] = f(x)
return memo[x]
return helper

@memoize
def is_prime(number):
# wheel / trial division
for i in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199):
if not (number % i):
return False
i = 203
while i <= number**0.5:
if not (number % i):
return False
i += 2
return True

digits_with_0 = string.digits
digits_without_0 = string.digits.replace("0", "")

def prime_neighbours(number):
n = str(number)
result = []
for i,v in enumerate(n):
if i == 0:
digits = digits_without_0
else:
digits = digits_with_0
for d in digits:
neighbour = list(n)
if d != neighbour[i]:
neighbour[i] = d
p = int("".join(neighbour))
if is_prime(p):
result.append(p)
return result

best = 100

def find_chain(current_chain, target, solutions):
global best
# if random.random() < 0.1:
# print(best, current_chain)
solutions = []
n = prime_neighbours(current_chain[-1])
# add hamming distance for sorting
z = [(sum(c1 == c2 for c1, c2 in zip(list(str(p1)), list(str(target)))), p1) for p1 in n]
z.sort(reverse=True)
for dummy, p2 in z:
if len(current_chain) + 1 < best and p2 not in current_chain:
if p2 == target:
solutions.append(current_chain + [p2])
print("#", len(solutions[-1]), solutions[-1])
best = len(current_chain)+1
else:
for s in find_chain(current_chain+[p2], target, solutions):
solutions.append(s)
return solutions

start = 71549
end = 10067

print(find_chain([start], end, []))
```