Enigmatic Code
Programming Enigma Puzzles
Enigma 45: Six squares – harder
Posted by on 1 December 2011
From New Scientist #1188, 3rd January 1980 [link] [link]
“You may remember,” said Mr Knull, “helping me recently to find the smallest possible set of three different integers P, Q and R such that P+Q, P+R, Q+R, P−Q, P−R and Q−R are all perfect squares. It was P=17, Q=8, R=−8.
“Here now is the harder version. Again you need to find P, Q and R satisfying the same conditions. But this time we don’t allow any of the three integers to be negative.”
Can you help me again, please, by finding the smallest set of three different integers which meet the conditions? By smallest I mean with P+Q+R as small as possible.
This puzzle refers back to Enigma 40.
A solution is given in New Scientist #1189, and in New Scientist #1190 it is noted that no correct solutions were recieved. However I think there is a solution where P+Q+R is smaller than that given as the solution in the magazine.
I sent New Scientist an email about this in December 2011, and in New Scientist #2853 an item was published in the Feedback section about my email, and I was awarded the prize for solving the puzzle 32 years late!
Note: This puzzle resurfaced as Project Euler Problem 142.
Note: This problem is known as Mengoli’s Six-Square Problem (1674). Euler found the smallest solution (1765), by hand.
[enigma45] [euler142]
Enigmatic Code 
It took me a couple of attempts before I found a satisfactory approach that would run in a sensible time. This Python program finds the published solution, and a “better” one in 1.6s (or 255ms using PyPy):
# consider: P+Q = a^2, P+R = b^2, Q+R = c^2 # a^2 + b^2 + c^2 = 2(P+Q+R) so we can minimise that # also: P > Q > R > 0 # so: a^2 > b^2 > c^2 # and also the following must be squares: # b^2 - c^2 = P-Q # a^2 - c^2 = P-R # a^2 - b^2 = Q-R from itertools import count from enigma import is_square, printf def solve(): for a in count(2): a2 = a ** 2 for b in range(1, a): b2 = b ** 2 if not is_square(a2 - b2): continue for c in range(0, b): c2 = c ** 2 if not is_square(a2 - c2): continue if not is_square(b2 - c2): continue (S, r) = divmod(a2 + b2 + c2, 2) if r > 0: continue (P, r) = divmod(a2 + b2 - c2, 2) if r > 0: continue (Q, r) = divmod(a2 + c2 - b2, 2) if r > 0: continue (R, r) = divmod(b2 + c2 - a2, 2) if r > 0: continue if not(P > Q > R > 0): continue yield (S, P, Q, R) # how many solutions to find import sys n = 2 if len(sys.argv) < 2 else eval(sys.argv[1]) for (S, P, Q, R) in solve(): printf("S={S} P={P} Q={Q} R={R}") if n == 1: break n -= 1Solution: P=434657, Q=420968, R=150568, P+Q+R=1006193.
This satisfies the conditions of the puzzle and is smaller than the published solution of: P=733025, Q=488000, R=418304, P+Q+R=1639329. (Which you can find by running with a command-line argument greater than 1).
Here’s my algorithm coded up in BBC BASIC, running on emulated hardware from the 80’s. For the details see [ http://jimpulse.blogspot.com/2012/02/retrocomputing-enigma-45.html ].
This slight modification makes the program run in half the time (800ms to find 2 solutions, or 213ms to find just the first one), by caching the (a², b²) pairs where a²−b² is a perfect square, and then using those to find the (b², c²) pairs. It also includes code to make sure the solutions come out in strict P+Q+R order (which slows it down again).
from itertools import count from collections import defaultdict from heapq import heappush, heappop from enigma import is_square, printf def solve(): # q: queue up the solutions q = [] # squares: a^2 -> list of b^2 where a^2 - b^2 is a perfect square squares = defaultdict(list) for a in count(2): a2 = a ** 2 for b in range(1, a): b2 = b ** 2 if not is_square(a2 - b2): continue squares[a2].append(b2) if not b2 in squares: continue for c2 in squares[b2]: if not is_square(a2 - c2): continue (S, r) = divmod(a2 + b2 + c2, 2) if r > 0: continue (P, r) = divmod(a2 + b2 - c2, 2) if r > 0: continue (Q, r) = divmod(a2 + c2 - b2, 2) if r > 0: continue (R, r) = divmod(b2 + c2 - a2, 2) if r > 0: continue if not(P > Q > R > 0): continue heappush(q, (S, P, Q, R)) # min S for this a ms = (a2 + 5) // 2 while q and q[0][0] < ms: yield heappop(q) # how many solutions to find import sys n = 2 if len(sys.argv) < 2 else eval(sys.argv[1]) for (S, P, Q, R) in solve(): printf("S={S} P={P} Q={Q} R={R}") if n == 1: break n -= 1I like programming and have just completed coding an algorithm of my own, to run in MATLAB. MATLAB is really far too powerful for the job! So how does an elapsed time of 13 minutes 38 seconds, for ‘a’ from 3 to 1000 sound?.
Will now look at the code above!
Tessa F
Here is an alternative Python implementation
from itertools import count, combinations from time import clock from sys import argv class Timer() : def __enter__(self): self.start = clock() def __exit__(self, *args): t = clock() - self.start if t < 10.0: print(' ({0:.2f} milliseconds)'.format(1000 * (clock() - self.start))) else: print(' ({0:.2f} seconds)'.format(clock() - self.start)) def is_square(x): t = int(x ** 0.5 + 0.1) return t if x == t * t else None def find(): # find Pythagorean triples where the hypotenuse (a) # can be represented in more than one way for a in count(2): sol = [] for x in range(1, a): t = a * a - x * x if t <= x * x: break y = is_square(t) if y: sol += [(y, x)] # take representations in pairs and allocate them to # (b, f) and (c, e) in such a way that b > c for cc in combinations(sol, 2): t1 = ( ((cc[0][0], cc[1][0]), (cc[1][1], cc[0][1])) if cc[0][0] > cc[1][0] else ((cc[1][0], cc[0][0]), (cc[0][1], cc[1][1])) ) t2 = ( ((cc[0][0], cc[1][1]), (cc[1][0], cc[0][1])) if cc[0][0] > cc[1][1] else ((cc[1][0], cc[0][1]), (cc[0][0], cc[1][1])) ) # check if these give a square for the final value (d) for (b, c), (e, f) in (t1, t2): d = is_square(b ** 2 - c ** 2) if d and len(set((a, b, c, d, e, f))) == 6: p, q, r = a * a + d * d, a * a - d * d, c * c - f * f if not ((p & 1) or (q & 1) or (r & 1)): yield p // 2, q // 2, r // 2, (p + q + r) // 2 n = 2 if len(argv) < 2 else eval(argv[1]) with Timer(): for p, q, r, s in find(): print('S ={0:8d} P ={1:8d} Q ={2:8d} R ={3:8d}'.format(s, p, q, r)) n -= 1 if not n: breakThanks for that. And for showing me the
[code]...[/code]tags – I didn’t know those worked in WordPress. Although with a bit of digging I’ve found you can use[sourcecode language="python"]...[/sourcecode]to get syntax highlighting too!An interesting find, Jim – can you edit the tags to add syntax highlighting to my post?
Yep. Done. I had a look to see if I could allow people to edit their own comments on wordpress.com, but I couldn’t find a setting for that.
% Another candidate for a MiniZinc solution % - but longer run-time than Enigma 40 include "globals.mzn"; var 0..500000: P; var 0..500000: Q; var 0..500000: R; predicate is_square(var int: c) = let { var 0..ceil(pow(int2float(ub(c)),(1/2.0))): z } in z * z = c ; constraint all_different([P,Q,R]); constraint is_square (P+Q) /\ is_square(P+R) /\ is_square(Q+R) /\ is_square(P-Q) /\ is_square(P-R) /\ is_square(Q-R); solve minimize(P+Q+R); output["P, Q, R = ",show(P)," , ",show(Q)," , ",show(R)]; % Running Enigma 45 Six Squares.mzn % P, Q, R = 434657 , 420968 , 150568 %---------- %Finished in 3m 46sTo identify the actual squares in this enigma:
As P, Q, R = 434657 , 420968 , 150568 :
1) P + Q = 855625 (925 * 925)
2) P + R = 585225 (765 * 765)
3) Q + R = 571536 (756 * 756)
4) P – Q = 13689 (117 * 117)
5) P – R = 284089 (533 * 533)
6) Q – R = 270400 (520 * 520)
That’s interesting. I did write a MiniZinc model to try and solve this some time ago, but I didn’t put upper bounds on the integers, and also used constraints to check for squares, like this:
When I gave it to the [[
mzn-gecode]] solver to chew on, it didn’t come back after 10 minutes, so I gave up on it.