## Euler’s Sum Of Powers Conjecture

### April 17, 2015

We’ll assume n = 5 and give three solutions of increasing efficiency. First the obvious brute-force solution:

```(define (euler n)
(do ((e 5 (+ e 1))) ((= e n))
(do ((d (- e 1) (- d 1))) ((zero? d))
(do ((c (- d 1) (- c 1))) ((zero? c))
(do ((b (- c 1) (- b 1))) ((zero? b))
(do ((a (- b 1) (- a 1))) ((zero? a))
(when (= (+ (expt a 5) (expt b 5)
(expt c 5) (expt d 5))
(expt e 5))
(display a) (display "^5 + ")
(display b) (display "^5 + ")
(display c) (display "^5 + ")
(display d) (display "^5 = ")
(display e) (display "^5")
(newline))))))))
> (time (euler 145))
27^5 + 84^5 + 110^5 + 133^5 = 144^5
(time (euler 145))
5919 collections
188949 ms elapsed cpu time, including 517 ms collecting
197025 ms elapsed real time, including 429 ms collecting
49846940832 bytes allocated, including 49849177312 bytes reclaimed
```

There are two sources of inefficiency in that program: recomputation of fifth-powers that have previously been calculated, and all the computations of all the loops. Our second version caches the fifth-power computations:

```(define (euler n)
(let ((cache (make-vector n)))
(do ((i 1 (+ i 1))) ((= i n))
(vector-set! cache i (expt i 5)))
(do ((e 5 (+ e 1))) ((= e n))
(do ((d (- e 1) (- d 1))) ((zero? d))
(do ((c (- d 1) (- c 1))) ((zero? c))
(do ((b (- c 1) (- b 1))) ((zero? b))
(do ((a (- b 1) (- a 1))) ((zero? a))
(when (= (+ (vector-ref cache a)
(vector-ref cache b)
(vector-ref cache c)
(vector-ref cache d))
(vector-ref cache e))
(display a) (display "^5 + ")
(display b) (display "^5 + ")
(display c) (display "^5 + ")
(display d) (display "^5 = ")
(display e) (display "^5")
(newline)))))))))
> (time (euler 145))
27^5 + 84^5 + 110^5 + 133^5 = 144^5
(time (euler 145))
5951 collections
104692 ms elapsed cpu time, including 264 ms collecting
108670 ms elapsed real time, including 407 ms collecting
50129778992 bytes allocated, including 50127622192 bytes reclaimed
```

An improvement of not quite two times. Our third implementation short-circuits the loops, stopping them when the requested sum is no longer possible:

```(define (euler n)
(let ((cache (make-vector n)))
(do ((i 1 (+ i 1))) ((= i n))
(vector-set! cache i (expt i 5)))
(do ((e 5 (+ e 1))) ((= e n))
(list-of (list a b c d e)
(d range (- e 1) 0)
(c range (- d 1) 0)
(c-sum is
(+ (vector-ref cache c)
(vector-ref cache d)))
(< c-sum (vector-ref cache e))
(b range (- c 1) 0)
(b-sum is
(+ c-sum
(vector-ref cache b)))
(< b-sum (vector-ref cache e))
(a range (- b 1) 0)
(= (+ b-sum (vector-ref cache a))
(vector-ref cache e))
(display a) (display "^5 + ")
(display b) (display "^5 + ")
(display c) (display "^5 + ")
(display d) (display "^5 = ")
(display e) (display "^5")
(newline)))))
(time (euler 145))
1298 collections
49421 ms elapsed cpu time, including 111 ms collecting
50556 ms elapsed real time, including 101 ms collecting
10934479232 bytes allocated, including 10933030352 bytes reclaimed
```

That’s about four times better than the original. But I’m not particularly happy with it; does anybody have any better ideas? Suggestions are welcome. You might also want to parameterize the size of n. You can run the program at http://ideone.com/RVcDBZ.

Pages: 1 2

### 7 Responses to “Euler’s Sum Of Powers Conjecture”

1. Paul said

In Python. I could not find anything clever.

```from time import clock
from itertools import combinations

def check(n, limit, first=2):
"""limitation: all terms must be different"""
pown = [i ** n for i in range(1, limit)]
pvals = set(pown)
for n_terms in range(first, n):
for p in combinations(pown, n_terms):
if sum(p) in pvals:
print p, clock() - start

def check5(limit):
pown = [i ** 5 for i in range(limit)]
pvals = set(pown)
for i, pi in enumerate(pown[1:], 1):
for j, pj in enumerate(pown[i:],i):
for k, pk in enumerate(pown[j:],j):
for l, pl in enumerate(pown[k:],k):
if pi+pj+pk+pl in pvals:
print i,j,k,l, clock() - start

start = clock()
check(5, 145) # time=4.85 secs
start = clock()
check5(145)   # time=3.44 secs
#Processor	Intel(R) Core(TM) i7-2630QM CPU @ 2.00GHz, 2001 MHz, 4 core
```
2. Paul said

A much faster method from the Lander and Parkin paper. The answer comes in a flash. I feel this can be improved, but this is already much faster.

```from __future__ import division

from time import clock
from bisect import bisect

def decompose(k, t, n):
pown = [i ** 5 for i in range(t)]
pvals = set(pown)
def find(k, t, n, data):
t5 = bisect(pown, t) - 1
tn5 = bisect(pown, t/n) - 1
if n == 2:
for v in range(tn5, t5+1):
if data and v > data[-1]:
continue
w = t - pown[v]
w5 = bisect(pown, w) - 1
if pown[w5] + pown[v] == t:
print "--", data, v, w5
else:
for s in range(t5, tn5 - 1, -1):
u = t - pown[s]
if data and data[-1] < s:
continue
find(k, u, n - 1, data + [s])
find(5, t ** 5, n, [])

start = clock()
for i in range(5, 145):
decompose(5, i, 4)
```
3. programmingpraxis said

@Paul: Please provide a citation to the paper.

4. Paul said

@ programmingpraxis: the paper: Lander and Parkin
I made a straightforward implementation. I feel it can be speeded up more. I did not use the bit at the bottom of page 101 with the mod 30 stuff.

5. Paul said

Another version in Python. It only produces “primitive” solutions, where the terms do not have a common factor. It produces the 2 tables from the Lander and Parkin paper. Also an example from Wikipedia for k=7 is found after a long time (about 40 minutes).

```from __future__ import division

from time import clock
from bisect import bisect_left, bisect_right
from fractions import gcd

def is_primitive(data):
"""if data is primitive, there is no common factor between all items"""
it = iter(data)
g = next(it, None)
if g:
for i in it:
g = gcd(g, i)
if g == 1:
return True
return False

def decompose_until(k, t, n, first=1):
""" decompose all numbers up to t as the sum of at most n terms of the form i ** k
"""
pown = [i ** k for i in range(t)]
pvals = set(pown)
L = len(pown)
def smin(target, n, hi=L):
"""returns the index in pown such that index ** k  > target // n
"""
t = target // n
return bisect_right(pown, t, hi=hi)
def smax(target, n, hi=L):
"""returns the index in pown such that index ** k  < target
if target is in pown then the index in pown is the returned index + 1
"""
return bisect_left(pown, target, hi=hi) - 1
def decompose(k, t, n):
"""decompose t**k as sum of terms of the form i ** k"""
tsafe = t
t = t ** k
Q = [(t - pown[s], n-1, [s]) for s in range(smin(t, n), smax(t, n)+1)]
while Q:
t, n, data = Q.pop()
smx = smax(t, n, data[-1]+1) + 1
if t in pvals and smx < data[-1] and is_primitive(data + [smx]):
assert sum(pown[i] for i in (data + [smx])) == tsafe ** k, "{:d} {:s}".format(tsafe, str(data+ [smx]))
print "**", tsafe, data + [smx], clock() - start
continue
smn = smin(t, n, data[-1]+1)
if n == 2:
for v in range(smn, smx + 1):
if v > data[-1]: continue
w = t - pown[v]
if w in pvals:
w5 = smax(w, n, data[-1]+1) + 1
if w5 <= v and is_primitive(data + [v, w5]):
assert sum(pown[i] for i in (data + [v, w5])) == tsafe ** k, "{:d} {:s}".format(tsafe, str(data+ [v, w5]))
print "--", tsafe, data + [v, w5], clock() - start
else:
for s in range(smx, smn - 1, -1):
if s > data[-1]: continue
Q.append((t - pown[s], n-1, data + [s]))
start = clock()
for i in range(first, t):
decompose(k, i, n)

decompose_until(5, 250, 4, first=1)
"""
k = 7 (Example wikipedia - takes long time - started @ 568
-- 568 [525, 439, 430, 413, 266, 258, 127] 2378.10256164 k=7
k = 5, n = 6 (Tabel 1 from Lander and Parkin)
-- 12 [11, 9, 7, 6, 5, 4] 0.361283899936
-- 30 [29, 19, 16, 11, 10, 5] 0.404110327133
-- 32 [28, 24, 22, 17, 16, 15] 0.430216805645
-- 67 [66, 34, 31, 29, 20, 7] 1.22579653624
-- 67 [66, 36, 31, 23, 18, 13] 1.24291212316
** 72 [67, 47, 46, 43, 19] 1.58158665811
-- 78 [64, 61, 58, 48, 35, 22] 2.25165772931
** 94 [84, 79, 37, 23, 21] 5.02510429763
-- 99 [96, 67, 20, 19, 13, 4] 6.27318537524
-- 99 [89, 73, 64, 60, 17, 6] 6.3922375882
k = 5, n = 5 (Tabel 2 from Lander and Parkin)
-- 72 [67, 47, 46, 43, 19] 0.211067094291
-- 94 [84, 79, 37, 23, 21] 0.604523412444
-- 107 [100, 80, 57, 43, 7] 1.00687568954
** 144 [133, 110, 84, 27] 3.26073688262
k = 5, n = 4
-- 144 [133, 110, 84, 27] 0.257228992944
"""
```
6. Marijn said

Somehow the superscript 5’s got left out of the 2 last terms, I found this in the html source: 133 = 144

7. programmingpraxis said

@Marijn: Fixed. Thanks.