Programming Enigma Puzzles
From New Scientist #1329, 28th October 1982 [link]
First, place the whole numbers from 1 to 16, one in each square. Next, multiply the numbers in each pair of touching squares and write the product in the little circle between the squares. Finally, add up the numbers in the circles. That gives you the “product-sum”, which you want to make as large as possible.
How big can you make it?
If you follow the Google Books link you will find that this Enigma was published next to a review (and advert) for the book Enigmas by Robert Eastaway, featuring 140 pages of Enigma puzzles selected from New Scientist. I have found a second hand copy of this and ordered it.
[enigma184]
 | Frits on Enigma 773: Duodecimal |
 | Jim Randell on Enigma 773: Duodecimal |
 | Brian Gladman on Enigma 773: Duodecimal |
| Frits on Enigma 773: Duodecimal |
| Frits on Enigma 773: Duodecimal |
| Frits on Enigma 773: Duodecimal |
| Frits on Enigma 773: Duodecimal |
| Jim Randell on Enigma 773: Duodecimal |
Enigmatic Code 
This is a directed search that starts from the observation that the numbers in the centre 4 squares are used in 4 products each and so are likely to be larger numbers, likewise the numbers in the 4 corner squares are only used in two products, so are likely to be smaller. (The remaining 8 edge squares are used in 3 products). It considers possible arrangements of the centre squares and calculates an upper bound for product sum that can be made using the remaining squares, if this is less than the best measure we’ve already found this arrangement is skipped.
It takes 28.2s to run to completion. Without the upper bound pruning I think it would take around 14 days to run to completion (although it does find the maximum solution in under 1s).
from itertools import combinations, permutations from enigma import irange, diff, printf # the numbers (largest first) ns0 = list(irange(16, 1, step=-1)) r = 0 # choose 4 numbers for the central 4 squares for m in combinations(ns0, 4): ns1 = diff(ns0, m) # let's put the largest number in position 10 and have the number at position 9 # greater than that at position 6 (to eliminate symmetrical solutions) for (n5, n6, n9, n10) in ((m[3], m[2], m[1], m[0]), (m[2], m[3], m[1], m[0]), (m[1], m[3], m[2], m[0])): # determine the central product s0 = (n5 * n6) + (n5 * n9) + (n6 * n10) + (n9 * n10) # choose small numbers for the corners for c in combinations(ns1[::-1], 4): # what's left are the edges ns2 = diff(ns1, c) # calculate an upper bound for the remaining products x = sorted(ns2) u = ((m[0] + c[3]) * (x[7] + x[6]) + (m[1] + c[2]) * (x[5] + x[4]) + (m[2] + c[1]) * (x[3] + x[2]) + (m[3] + c[0]) * (x[1] + x[0]) + ((x[7] * x[6]) + (x[5] * x[4]) + (x[3] * x[2]) + (x[1] * x[0]))) # if it's less then our current best skip it if s0 + u < r: continue # permute the corners for (n0, n3, n12, n15) in permutations(c): # what's left are the edges for e in permutations(ns2): (n1, n2, n4, n7, n8, n11, n13, n14) = e # compute the product sum s = s0 + ((n0 * n1) + (n1 * n2) + (n2 * n3) + (n4 * n5) + (n6 * n7) + (n8 * n9) + (n10 * n11) + (n12 * n13) + (n13 * n14) + (n14 * n15) + (n0 * n4) + (n4 * n8) + (n8 * n12) + (n1 * n5) + (n9 * n13) + (n2 * n6) + (n10 * n14) + (n3 * n7) + (n7 * n11) + (n11 * n15)) if s > r: r = s printf("[{r}] {n0} {n1} {n2} {n3} / {n4} {n5} {n6} {n7} / {n8} {n9} {n10} {n11} / {n12} {n13} {n14} {n15}") printf("max measure = {r}")Solution: The maximum possible product sum is 2345.
Here’s a diagram of the arrangement:
Labelling the squares a0,…a15, across then down, and colouring them as in a chessboard, if all the black squares, a0,a2,a5,a7,a8,a10,a13,a15 are populated, we can deduce the maximum product-sum possible by populating the white squares with the remaining numbers by using the rearrangement inequality.
This program considers permutations of numbers populating the black squares, eliminating reflections, and eliminating permutations where the black squares total more than the white squares. It still runs slower than Jim’s version, about 3 minutes under PyPy.
from itertools import permutations as perm, combinations as comb, starmap,izip from operator import mul best=0 set1to16=frozenset(range(1,17)) for a15,a0 in comb(range(1,17),2): for a8,a2 in comb(set1to16-set((a15,a0)),2): for a5,a7,a10,a13 in perm(set1to16-set((a15,a0,a2,a8)),4): if a0+a2+a5+a7+a8+a10+a13+a15 <= 68: m = sorted((a0+a5+a2,a2+a7,a0+a5+a8,a2+a5+a7+a10,a5+a8+a10+a13,a10+a7+a15,a8+a13,a10+a13+a15)) r = sorted(set1to16-set((a0,a2,a5,a7,a8,a10,a13,a15))) score = sum(starmap(mul,izip(m,r))) # dot product of m and r if score>best: print "best score",score best=scoreSome optimisation gets the run time down to 46 sec
from itertools import permutations as perm, combinations as comb, starmap,izip from operator import mul best=0 set1to16=frozenset(range(1,17)) for a15,a0 in comb(range(1,17),2): for a8,a2 in comb(set1to16-set((a15,a0)),2): for x in comb(set1to16-set((a15,a0,a2,a8)),4): if a0+a15+a8+a2 + sum(x) <= 68: sx=sorted(x) m=sorted((a0+sx[3]+a2,a2+sx[3],a0+sx[3]+a8,a2+sx[3]+sx[2]+sx[1],sx[3]+a8+sx[2]+sx[1],sx[3]+sx[2]+a15,a8+sx[3],sx[3]+sx[2]+a15)) r = sorted(set1to16-set((a0,a15,a8,a2))-set(x)) upperbound = sum(starmap(mul,izip(m,r))) # dot product of m and r if upperbound > best: for a5,a7,a10,a13 in perm(x): m = sorted((a0+a5+a2,a2+a7,a0+a5+a8,a2+a5+a7+a10,a5+a8+a10+a13,a10+a7+a15,a8+a13,a10+a13+a15)) score = sum(starmap(mul,izip(m,r))) # dot product of m and r if score>best: print "best score",score best=score