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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: