In a blog post, John Baez proposes the following puzzle:

Consider solutions of

with positive integers and show that the largest possible value of is .

We can start by writing the following simple Python program to solve the problem, using integer arithmetic to avoid floating point problems:

from itertools import count for r in count(1): for q in xrange(1, r + 1): for p in xrange(1, q + 1): # 1/p + 1/q + 1/r = 1/2 # qr + pr + pq = pqr/2 # 2(qr + pr + pq) = pqr if 2*(q*r + p*r + p*q) == p*q*r: print p, q, r

This program shows us ten possible values (including 42 in the last triple),

6 6 6 4 8 8 5 5 10 4 6 12 3 12 12 3 10 15 3 9 18 4 5 20 3 8 24 3 7 42

but, as it’s trying to enumerate all integer triples, it doesn’t halt.

To be sure not to miss any interesting integer triple, we must explore them more carefully:

from itertools import count for p in count(1): # checks if we are leaving room for 1/q + 1/r # 1/p >= 1/2 if 2 >= p: continue # we need smaller values of 1/p # checks if we can reach 1/2 with the biggest legal values for 1/q and 1/r # 1/p + 1/p + 1/p < 1/2 if 3 * 2 < p: break for q in count(p): # checks if we are leaving room for 1/r # 1/p + 1/q >= 1/2 if 2 * (q + p) >= p * q: continue # checks if we can reach 1/2 with the biggest legal value for 1/r # 1/p + 1/q + 1/q < 1/2 if 2 * (q + 2 * p) < p * q: break for r in count(q): lhs = 2 * (q * r + p * r + p * q) rhs = p * q * r # 1/p + 1/q + 1/r > 1/2 if lhs > rhs: continue # 1/p + 1/q + 1/r < 1/2 elif lhs < rhs: break # 1/p + 1/q + 1/r = 1/2 else: print p, q, r

The results are the same, but now we can be sure of having the whole solution.