Programming Praxis

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

1. Paul said

   April 17, 2015 at 2:21 PM

   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

   April 17, 2015 at 3:22 PM

   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

   April 17, 2015 at 6:48 PM

   @Paul: Please provide a citation to the paper.
4. Paul said

   April 17, 2015 at 7:53 PM

   @ 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

   April 18, 2015 at 1:45 PM

   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

   May 8, 2015 at 7:33 AM

   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

   May 8, 2015 at 6:41 PM

   @Marijn: Fixed. Thanks.