Notebook

Euler's Sum of Powers Conjecture¶

In 1769, Leonhard Euler conjectured that for all integers n and k greater than 1, if the sum of n kth powers of positive integers is itself a kth power, then n is greater than or equal to k. For example, this would mean that no sum of a pair of cubes (a³ + b³) can be equal to another cube (c³), but a sum of three cubes can, as in 3³ + 4³ + 5³ = 6³.

It took 200 years to disprove the conjecture: in 1966 L. J. Lander and T. R. Parkin published a refreshingly short article giving a counterexample of four fifth-powers that summed to another fifth power. They found it via a program that did an exhaustive search. Can we duplicate their work and find integers greater than 1 such that a⁵ + b⁵ + c⁵ + d⁵ = e⁵ ?

Algorithm¶

An exhaustive O(m⁴) algorithm woud be to look at all values of a, b, c, d < m and check if a⁵ + b⁵ + c⁵ + d⁵ is a fifth power. But we can do better: a sum of four numbers is a sum of two pairs of numbers, so we are looking for

pair₁ + pair₂ = e⁵ where pair₁ = a⁵ + b⁵ and pair₂ = c⁵ + d⁵.

We will define pairs be a dict of {a⁵ + b⁵: (a⁵, b⁵)} entries for all a ≤ b < m; for example, for a=2 and b=10, the entry is {100032: (32, 100000)}. Then we can ask for each pair₁, and for each e, whether there is a pair₂ in the dict that makes the equation work. There are O(m²) pairs and O(m) values of e, and dict lookup is O(1), so the whole algorithm is O(m³):

In [1]:

import itertools

def euler(m):
    """Yield tuples (a, b, c, d, e) such that a^5 + b^5 + c^5 + d^5 = e^5,
    where all are integers, and 1 < a ≤ b ≤ c ≤ d < e < m."""
    powers = [e**5 for e in range(2, m)] 
    pairs  = {sum(pair): pair 
              for pair in itertools.combinations_with_replacement(powers, 2)}
    for pair1 in pairs:
        for e5 in powers:
            pair2 = e5 - pair1
            if pair2 in pairs:
                yield fifthroots(pairs[pair1] + pairs[pair2] + (e5,))
    
def fifthroots(nums): 
    "Sorted integer fifth roots of a collection of numbers." 
    return tuple(sorted(int(round(x ** (1/5))) for x in nums))

Let's look for a solution (arbitrarily choosing m=500):

In [2]:

%time next(euler(500))

CPU times: user 1.07 s, sys: 21.4 ms, total: 1.09 s
Wall time: 1.11 s

Out[2]:

(27, 84, 110, 133, 144)

That was easy, and it turns out this is the same answer that Lander and Parkin got: 27⁵ + 84⁵ + 110⁵ + 133⁵ = 144⁵.

We can keep going, collecting all the solutions up to *m*=1000:

In [3]:

%time set(euler(1000))

CPU times: user 1min 53s, sys: 706 ms, total: 1min 54s
Wall time: 1min 57s

Out[3]:

{(27, 84, 110, 133, 144),
 (54, 168, 220, 266, 288),
 (81, 252, 330, 399, 432),
 (108, 336, 440, 532, 576),
 (135, 420, 550, 665, 720),
 (162, 504, 660, 798, 864)}

All the answers are multiples of the first one (this is easiest to see in the middle column: 110, 220, 330, ...). Since 1966 other mathematicians have found other solutions, but all we need is one to disprove Euler's conjecture.