Teaser 3177: Prime birthday card
From The Sunday Times, 13th August 2023 [link] [link]
On her 11th birthday, Ann’s older brother Ben gave her a card on which he had drawn a 5×5 square grid with a different prime number in each of the cells. Every 5-cell row, column and diagonal had the same sum and, noting that 1 is not a prime, he used primes that made this sum as small as possible.
The centre cell contained Ann’s age. All but one of the primes containing a digit 7 were on the two diagonals, the largest of these being in the bottom right corner. Two cells in the far right column contained a digit 5 as did two cells of the 4th row down. Four of the cells in the middle row contained a digit 3, and the largest prime on the grid was in the same column as two single digit primes.
What are the five prime numbers on the top row of the card?
[teaser3177]
Hugo, 
NigelR, and 1 other are discussing. 
Frits, 
GeoffR, and 1 other are discussing. 
S2T2 
Jim Randell 5:42 pm on 11 August 2023 Permalink |
We can exclude 2 from the list of primes, as we need the sum of all 12 magic lines to have the same parity (as they are all equal to the magic constant).
The magic constant of the square is the sum of the numbers used, divided by 5, so we can look for candidate sets of primes that give the smallest magic constant.
There are two potential sets of 25 primes that would give a magic constant of 233:
And it is possible to construct a Magic Square using either of these sets, so one of these sets must be the set that Ben used.
(I used a separate program to find the minimal sets, see: [@replit]).
In the completed grid the centre cell is occupied by 11, so there are 8 remaining cells on the diagonals.
The first of these sets has 8 numbers that include the digit 7. So is possible for 7 of the remaining diagonal cells to use 7 of these numbers, and the remaining number is in one of the non-diagonal cells.
The second of these sets has 10 numbers that include the digit 7. So, even if all 8 remaining diagonal cells are filled with one of these numbers there would be 2 numbers left over, so this set is not viable.
So the set of primes used by Ben must be set (1).
We can now look for ways to fill out the grid as specified by the remaining requirements.
I used a MiniZinc program to generate candidate magic squares that fulfil some (most) of the requirements, and then perform additional checks in Python to find arrangements that satisfy all the requirements. I used the minizinc.py library for this.
This program performs an exhaustive search in 702ms.
from enigma import (primes, ediv, nsplit, diff, chunk, unzip, implies, icount, lt, join, int2base, printf) from minizinc import (MiniZinc, var) # numbers to use (= set (1)) numbers = diff(primes.between(3, 103), {97}) # the magic constant k = ediv(sum(numbers), 5) # find numbers that include the digit 3, 5, 7 (n3s, n5s, n7s) = (set(), set(), set()) for n in numbers: ds = nsplit(n) if 3 in ds: n3s.add(n) if 5 in ds: n5s.add(n) if 7 in ds: n7s.add(n) # format a sequence for MiniZinc fseq = lambda xs, sep=", ", enc="{}": join(xs, sep=sep, enc=enc) # construct the model: # # A B C D E # F G H I J # K L M N P # Q R S T U # V W X Y Z # vs = "ABCDEFGHIJKLMNPQRSTUVWXYZ" model = f""" include "globals.mzn"; % possible numbers set of int: Number = {fseq(numbers)}; % decision variables for the cells of the square {var("Number", vs)}; % general constraints: % each cell is different constraint all_different({fseq(vs, enc="[]")}); % magic lines have the same sum (= k) constraint all_equal([{k}, A+B+C+D+E, F+G+H+I+J, K+L+M+N+P, Q+R+S+T+U, V+W+X+Y+Z, % rows A+F+K+Q+V, B+G+L+R+W, C+H+M+S+X, D+I+N+T+Y, E+J+P+U+Z, % cols A+G+M+T+Z, E+I+M+R+V % diags ]); % some specific constraints: % centre cell is 11 constraint M = 11; % and the other 4 cells in the middle row (K, L, N, P) contain a digit 3 constraint {{K, L, N, P}} subset {fseq(n3s)}; % 2 cells in the far right column contain a digit 5 constraint card({{E, J, P, U, Z}} intersect {fseq(n5s)}) = 2; % as did 2 cells in the 4th row down constraint card({{Q, R, S, T, U}} intersect {fseq(n5s)}) = 2; % all but one of the primes containing a '7' are on the diagonals constraint card({fseq(n7s)} diff {{A, E, G, I, R, T, V, Z}}) = 1; % and the largest of these is Z constraint max({fseq(n7s)} intersect {{A, E, G, I, R, T, V, Z}}) = Z; solve satisfy; """ # load the model p = MiniZinc(model, solver="minizinc --solver gecode --all", verbose=0) # and solve it for s in p.solve(): ns = tuple(s[k] for k in vs) rows = list(chunk(ns, 5)) cols = list(unzip(rows)) (A, B, C, D, E, F, G, H, I, J, K, L, M, N, P, Q, R, S, T, U, V, W, X, Y, Z) = ns # centre cell is 11 if not (M == 11): continue # and the other 4 cells in the middle row contain a digit 3 if not n3s.issuperset([K, L, N, P]): continue # 2 cells in the far right column have a digit 5 if not (len(n5s.intersection([E, J, P, U, Z])) == 2): continue # also 2 of the cells in the 4th row down if not (len(n5s.intersection([Q, R, S, T, U])) == 2): continue # all but one of the numbers containing the digit 7 are on the diagonals x7s = diff(n7s, {A, E, G, I, R, T, V, Z}) if not (len(x7s) == 1): continue # and the largest of these numbers is Z if not (max(diff(n7s, x7s)) == Z): continue # the largest prime on the grid is in the same column as 2 single digit primes maxn = max(ns) if not all(implies(maxn in col, icount(col, lt(10)) == 2) for col in cols): continue # output the grid fmt = lambda n: int2base(n, width=3, pad=' ') for r in rows: printf("{r}", r=join(r, sep=" ", enc="[]", fn=fmt)) printf()Solution: The numbers on the top row of the card are: 17, 7, 101, 41, 67.
The complete grid is:
LikeLike
Frits 9:11 am on 12 August 2023 Permalink |
Thanks. I wondered how you would do the y contains ‘x’ thing in MiniZinc without string support.
LikeLike
Jim Randell 10:55 am on 12 August 2023 Permalink |
Originally I had fewer constraints in the MiniZinc model (which is why they are all tested in Python), but it ran faster if the constraints were implemented in the model (apart from the “largest prime” constraint, which was messier to do in MiniZinc, and also slowed down the overall execution).
As it is the model now only generates a single candidate, so the cost of verifying all the requirements in Python is negligible. (Although the repeated tests could be removed to give a shorter program).
But I used a restricted version of the model to find the minimal sets of primes that form a Magic Square.
LikeLike
Frits 10:09 am on 13 August 2023 Permalink |
from itertools import combinations, permutations # decompose <t> into <k> numbers from <ns> def decompose(t, k, ns, s=tuple()): if k == 1: if t in ns: yield s + (t, ) else: for n in ns: if t - n < minPr: continue yield from decompose(t - n, k - 1, ns.difference([n]), s + (n ,)) # decompose <t> into <k> increasing odd prime numbers from <ns> def decompose_inc(t, k, ns, s=[]): if k == 1: if t in ns: yield set(s + [t]) else: for i, n in enumerate(ns): if t < k * n + k * (k - 1): break yield from decompose_inc(t - n, k - 1, ns[i + 1:], s + [n]) # exclude entries that occur in <s> ps = lambda s: {p for p in Pr if p not in s} # A B C D E # F G H I J # K L M N O # P Q R S T # U V W X Y def solve(): sol = 0 # all but one of the primes containing a digit 7 were on two diagonals # choose <no7s - 1> 7's for the diagonals for c7 in combinations(sevens, no7s - 1): rest = 2 * s - 2 * M - sum(c7) # rest sum for other diagonals Y = c7[-1] # largest of diagonal 7's being in the bottom right corner match(no7s): case 9: oth = ((), ) # tuple of an empty tuple case 8: oth = ((rest, ), ) if rest in Pr and '7' not in str(rest) \ else tuple() case other: oth = tuple(decompose(rest, 9 - no7s, {p for p in Pr if '7' not in str(p)})) if not oth: continue # for all combinations of diagonal cells without '7' (or center) for o in oth: # permute unknown diagonal values for A, G, S, E, I, Q, U in permutations(c7[:-1] + o): # only check one diagonal sum if A + G + M + S + Y != s: continue s1 = {A, G, M, S, Y, E, I, Q, U} # four of the cells in the middle row contained a digit 3 for K, L, N, O in decompose(s - M, 4, {p for p in Pr.difference(s1) if '3' in str(p)}): # determining J, T followed by determining P, R is a lot faster # than determining P, R, T and calculating J (credit: B. Gladman) # 5th column for J, T in decompose(s - E - O - Y, 2, ps(s2 := s1 | {K, L, N, O})): # two cells in the far right column contained a digit 5 if sum('5' in str(x) for x in [E, J, O, T, Y]) != 2: continue # 4th row for P, R in decompose(s - Q - S - T, 2, ps(s3 := s2 | {J, T})): # two cells in the 4th row down contained a digit 5 if sum('5' in str(x) for x in [P, Q, R, S, T]) != 2: continue F = s - (A + K + P + U) # 1st column if F in s3 | {P, R} or F not in Pr: continue H = s - (F + G + I + J) # 2nd row if H in (s4 := s3 | {P, R, F}) or H not in Pr: continue # 2nd column for B, V in decompose(s - G - L - Q, 2, ps(s5 := s4 | {H})): # 3rd column for C, W in decompose(s - H - M - R, 2, ps(s6 := s5 | {B, V})): D = s - (A + B + C + E) # 1st row X = s - (U + V + W + Y) # 5th row if any(x not in Pr for x in [D, X]): continue if len(s6 | {C, W, D, X}) != 25: continue mat = [(A, B, C, D, E), (F, G, H, I, J), (K, L, M, N, O), \ (P, Q, R, S, T), (U, V, W, X, Y)] # transpose matrix tmat = list(zip(*mat)) # largest prime on the grid was in the same column as two # single digit primes if sum(max(Pr) in x for x in tmat \ if sum(y < 10 for y in x) >= 2) != 1: continue # double check sums for rows and columns if any(sum(x) != s for m in [mat, tmat] for x in m): continue print() for m in mat: print(''.join([f"{str(x):>4}" for x in m])) print("\nanswer:", mat[0]) # break if we have finished processing this magical sum <s> sol = 1 return sol # list of prime numbers up to 997 primes = [3, 5, 7] primes += [x for x in range(11, 100, 2) if all(x % p for p in primes)] primes += [x for x in range(101, 1000, 2) if all(x % p for p in primes)] total = sum(primes[:25]) # find first entry higher than total for which total = s * 5 for an odd s while total % 5 or total % 2 == 0: total += 1 M = 11 # the centre cell contained Ann's age primes = [p for p in primes if p != M] sol = 0 # potentially process all odd magical sums <s> up to a high number for s in range(total // 5, 10000, 2): # choose 24 prime numbers that (together with <M>) sum up to 5 * s for Pr in decompose_inc(5 * s - M, 24, primes): sevens = [p for p in Pr if '7' in str(p)] # 9 diagonal cells from which M is 11 if (no7s := len(sevens)) > 9: continue minPr = min(Pr) # solve the puzzle sol += solve() if sol: breakLikeLike
Frits 10:10 am on 13 August 2023 Permalink |
from enigma import SubstitutedExpression # decompose <t> into <k> increasing numbers from <ns> def decompose_inc(t, k, ns, s=[]): if k == 1: if t in ns: yield set(s + [t]) else: for i, n in enumerate(ns): if t < k * n + 2: break yield from decompose_inc(t - n, k - 1, ns[i + 1:], s + [n]) ''' A B C D E F G H I J K L M N O P Q R S T U V W X Y ''' # list of prime numbers primes = [3, 5, 7] primes += [x for x in range(11, 100, 2) if all(x % p for p in primes)] primes += [x for x in range(101, 1000, 2) if all(x % p for p in primes)] total = sum(primes[:25]) # find first entry higher than total for which total = s * 5 for an odd s while total % 5 or total % 2 == 0: total += 1 sol = 0 # potentially process all odd magical sums <s> up to a high number for s in range(total // 5, 10000, 2): for Pr in decompose_inc(5 * s, 25, primes): sevens = [p for p in Pr if '7' in str(p)] # 9 diagonal cells from which M is 11 if (no7s := len(sevens)) > 9: continue exprs = [] # diagonal exprs += ["A + G + M + S + Y == " + str(s)] exprs += ["sum('7' in str(x) for x in [A, G, M, S]) >= 2"] exprs += ["max(x for x in [A, G, M, S, Y] if '7' in str(x)) == Y"] exprs += ["E + I + M + Q + U == " + str(s)] # check for <no7s - 1> 7's on the diagonals exprs += ["sum('7' in str(x) for x in [A, G, M, S, Y, E, I, Q, U]) == " +\ str(no7s - 1)] exprs += ["max(x for x in [A, G, M, S, Y, E, I, Q, U] \ if '7' in str(x)) == Y"] # 3rd row exprs += ["K + L + M + N + O == " + str(s)] # four of the cells in the middle row contained a digit 3 exprs += ["sum('3' in str(x) for x in [K, L, M, N, O]) == 4"] # 5th column exprs += ["E + J + O + T + Y == " + str(s)] # two cells in the far right column contained a digit 5 exprs += ["sum('5' in str(x) for x in [E, J, O, T, Y]) == 2"] # 4th row exprs += ["P + Q + S + T + R == " + str(s)] # two cells in the 4th row down contained a digit 5 exprs += ["sum('5' in str(x) for x in [P, Q, R, S, T]) == 2"] # remaining horizontal exprs += ["A + B + C + D + E == " + str(s)] exprs += ["F + G + H + I + J == " + str(s)] exprs += ["U + V + W + X + Y == " + str(s)] # remaining vertical exprs += ["A + F + K + U + P == " + str(s)] exprs += ["B + G + L + Q + V == " + str(s)] exprs += ["C + H + M + R + W == " + str(s)] exprs += ["D + I + N + S + X == " + str(s)] # largest prime on the grid was in the same column # as two single digit primes exprs += ["sum(" + str(max(Pr)) + " in x for x in [(A, F, K, P, U), \ (B, G, L, Q, V), (C, H, M, R, W), (D, I, N, S, X), \ (E, J, O, T, Y)] if sum(y < 10 for y in x) >= 2) == 1"] # the alphametic puzzle p = SubstitutedExpression( exprs, answer="((A, B, C, D, E), (F, G, H, I, J), (K, L, M, N, O), \ (P, Q, R, S, T), (U, V, W, X, Y))", s2d=dict(M=11), base=max(Pr) + 1, digits=Pr, #denest=8, # use denest if not run with PyPy verbose=0, # use 256 to see the generated code ) # print answers for ans in p.answers(): sol = 1 print() for a in ans: print(''.join([f"{str(x):>4}" for x in a])) print(f"\nanswer: {ans[0]}") if sol: breakLikeLike
NigelR 11:28 am on 13 August 2023 Permalink |
Hi Frits
I could be missing something, and it may not be relevant to this problem, but the ‘break’ condition (line 10) in your function decompose_inc seems to risk missing some valid decompositions. For example, if you call decompose_inc(10,3,[1,2,3,4,5]), it doesn’t yield {1,4,5}. If you change the condition to ‘if t < k * n + 1: break' it seems to work correctly. Is there a reason for using +2 rather than +1?
LikeLike
Frits 2:10 pm on 13 August 2023 Permalink |
Hi Nigel, I chose the “+ 2” as there always is a difference of at least 2 between odd prime numbers.
LikeLike
NigelR 5:46 pm on 13 August 2023 Permalink |
Thanks Frits. I can now see that it is tailored to the context rather than intended as a generic decompose generator, which is how I was trying to use it!
LikeLike
Frits 2:09 pm on 13 August 2023 Permalink |
A MiniZinc solution with a matrix based on Jim’s program and Brian Gladman’s initial setup.
It runs in 411ms.
from enigma import (nsplit, join) from minizinc import MiniZinc # decompose <t> into <k> increasing numbers from <ns> def decompose_inc(t, k, ns, s=[]): if k == 1: if t in ns: yield set(s + [t]) else: for i, n in enumerate(ns): if t < k * n + 2: break yield from decompose_inc(t - n, k - 1, ns[i + 1:], s + [n]) # format a sequence for MiniZinc fseq = lambda xs, sep=", ", enc="{}": join(xs, sep=sep, enc=enc) # list of prime numbers up to 997 primes = [3, 5, 7] primes += [x for x in range(11, 100, 2) if all(x % p for p in primes)] primes += [x for x in range(101, 1000, 2) if all(x % p for p in primes)] total = sum(primes[:25]) # find first entry higher than total for which total = s * 5 for an odd s while total % 5 or total % 2 == 0: total += 1 sol = 0 # potentially process all odd magical sums <s> up to a high number for s in range(total // 5, 10000, 2): # choose 25 prime numbers that sum up to 5 * s for Pr in decompose_inc(5 * s, 25, primes): sevens = [p for p in Pr if '7' in str(p)] # 9 diagonal cells from which center is 11 if (no7s := len(sevens)) > 9: continue # find numbers that include the digit 3, 5, 7 (n3s, n5s, n7s) = (set(), set(), set()) for n in Pr: ds = nsplit(n) if 3 in ds: n3s.add(n) if 5 in ds: n5s.add(n) if 7 in ds: n7s.add(n) magic = s model = f""" include "globals.mzn"; % possible numbers set of int: Number = {fseq(Pr)}; set of int: SIZE = 1..5; array[SIZE, SIZE] of var Number: sq; %var set of int: r3 = array2set(row(sq, 3)); %var set of int: r4 = array2set(row(sq, 4)); %var set of int: c5 = array2set(col(sq, 5)); % the centre cell contained Ann's age constraint sq[3, 3] = 11; constraint alldifferent(i in SIZE, j in SIZE) (sq[i, j]); % same row sums constraint forall(i in SIZE) (sum(j in SIZE) (sq[i, j]) = {magic}); % same column sums constraint forall(j in SIZE) (sum(i in SIZE) (sq[i, j]) = {magic}); % same diagonal sums constraint sum(i in SIZE) (sq[i, i]) = {magic}; constraint sum(i in SIZE) (sq[i, 6 - i]) = {magic}; % all primes but one with '7' are located on the diagonals % and the largest of these is sq[5, 5] constraint sum(i, j in SIZE where i == j \/ i + j == 6) (sq[i, j] in {fseq(n7s)} /\ sq[i, j] <= sq[5, 5]) == {no7s - 1}; % and the other 4 cells in the middle row contain a digit 3 % constraint card(r3 intersect {fseq(n3s)}) = 4; constraint sum(j in SIZE) (sq[3, j] in {fseq(n3s)}) = 4; % 2 cells in the far right column contain a digit 5 %constraint card(c5 intersect {fseq(n5s)}) = 2; constraint sum(i in SIZE) (sq[i, 5] in {fseq(n5s)}) = 2; % as did 2 cells in the 4th row down %constraint card(r4 intersect {fseq(n5s)}) = 2; constraint sum(j in SIZE) (sq[4, j] in {fseq(n5s)}) = 2; solve satisfy; """ # load the model p = MiniZinc(model, solver="minizinc --solver gecode --all", verbose=0) # and solve it for s in p.solve(): sq = s['sq'] for y in sq: print(y) print("\nanswer:", sq[0]) sol = 1 if sol: breakLikeLike
Frits 5:03 pm on 13 August 2023 Permalink |
The largest prime requirement is not needed to find a unique solution.
Here is the MiniZinc constraint (it is faster to check in afterwards in Python).
% largest prime on the grid was in the same column as two single digit primes constraint sum(j in SIZE) ({max(Pr)} in col(sq, j) /\ card(array2set(col(sq, j)) intersect {{3, 5, 7}}) >= 2) == 1;LikeLike
Hugo 11:19 am on 20 August 2023 Permalink |
OEIS no. A164843 confirms that 233 is the smallest sum for an order-5 prime magic square.
Is it possible to arrange the same twenty-five primes (alternatively with 97 and 107 instead of 101 and 103) to form a ‘panmagic’ square, i.e where the broken diagonals also sum to 233?
LikeLike
Jim Randell 5:50 pm on 20 August 2023 Permalink |
I added these constraints into my MiniZinc program, and it has been running some for some time to see if it is possible to find a construction (which makes me suspect it is not). But I will keep you posted if it finds anything.
LikeLike
Hugo 6:38 pm on 20 August 2023 Permalink |
Thanks, Jim.
I too suspect it won’t work with sum 233. On the tangled web I found one with sum 395. There may be one with a smaller sum than that, but my program would take about 1000 years to find it.
LikeLike