Euler's sum of powers conjecture
You are encouraged to solve this task according to the task description, using any language you may know.
There is a conjecture in mathematics that held for over 200 years before it was disproved by the finding of a counterexample in 1966 by Lander and Parkin.
- Euler's (disproved) sum of powers conjecture
- At least k positive kth powers are required to sum to a kth power, except for the trivial case of one kth power: yk = yk.
Lander and Parkin are known to have used a brute-force search on a CDC 6600 computer restricting numbers to those less than 250.
- Task
Write a program to search for an integer solution to:
x05 + x15 + x25 + x35 == y5
Where all xi's and y are distinct integers between 0 and 250 (exclusive).
Show an answer here.
C
The trick to speed up was the observation that for any x we have x^5=x modulo 2, 3, and 5, according to the Fermat's little theorem. Thus, based on the Chinese Remainder Theorem we have x^5==x modulo 30 for any x. Therefore, when we have computed the left sum s=a^5+b^5+c^5+d^5, then we know that the right side e^5 must be such that s==e modulo 30. Thus, we do not have to consider all values of e, but only values in the form e=es+30k, for some starting value es, and any k. Also, we follow the constraints 1<=a<b<c<d<e<N in the main loop. <lang c>#include <stdio.h>
- include <time.h>
typedef long long mylong;
void compute(void) { const int N = 250; /* variables are from [1..N-1] */ const int M = 30; /* x^5 == x modulo M=2*3*5 */ int a, b, c, d, e, es; mylong p5[N+M], s, t, max, i;
for(i=0; i<N; ++i) p5[i] = i*i, p5[i]*=p5[i]*i; for(max = p5[N-1]; i<(N+M); p5[i++]=max+1);
for(a=1; a<N; ++a) for(b=a+1; b<N; ++b) for(c=b+1, i=p5[a]+p5[b]; c<N; ++c) for(d=c+1, es=d+((t=p5[c]+i)%M); ((s = t + p5[d])<=max) && p5[es]<=s; ++d) { for(e=es++; p5[e]<s; e+=M); /* jump over M=30 values for e>d */ if(p5[e]==s) { printf("%d %d %d %d %d\r\n", a, b, c, d, e); return; } } }
int main(void) { int tm = clock(); compute(); printf("time=%d milliseconds\r\n", (int)((clock()-tm)*1000/CLOCKS_PER_SEC)); return 0; }</lang>
- Output:
27 84 110 133 144 time=70 milliseconds
C++
<lang cpp>#include <iostream>
- include <cmath>
- include <set>
using namespace std;
int main() { const auto MAX = 250; double pow5[MAX]; set<double> pow5s; for (auto i = 1; i < MAX; i++) { pow5[i] = (double)i * i * i * i * i; pow5s.insert(pow5[i]);} for (auto x0 = 1; x0 < MAX; x0++) { for (auto x1 = 1; x1 < x0; x1++) { for (auto x2 = 1; x2 < x1; x2++) { for (auto x3 = 1; x3 < x2; x3++) { auto sum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3]; if (pow5s.find(sum) != pow5s.end()) { cout << x0 << " " << x1 << " " << x2 << " " << x3 << " " << pow(sum, 1.0 / 5.0) << endl; goto exit;}}}}}
exit:
return 0;} </lang>
- Output:
133 110 84 27 144
EchoLisp
To speed up things, we search for x0, x1, x2 such as x0^5 + x1^5 + x2^5 = a difference of 5-th powers. <lang lisp> (define dim 250)
- speed up n^5
(define (p5 n) (* n n n n n)) (remember 'p5) ;; memoize
- build vector of all y^5 - x^5 diffs - length 30877
(define all-y^5-x^5 (for*/vector [(x (in-range 1 dim)) (y (in-range (1+ x) dim))] (- (p5 y) (p5 x))))
- sort to use vector-search
(begin (vector-sort! < all-y^5-x^5) 'sorted)
;; find couple (x y) from y^5 - x^5
(define (x-y y^5-x^5) (for*/fold (x-y null) [(x (in-range 1 dim)) (y (in-range (1+ x ) dim))] (when (= (- (p5 y) (p5 x)) y^5-x^5) (set! x-y (list x y)) (break #t)))) ; stop on first
- search
(for*/fold (sol null) [(x0 (in-range 1 dim)) (x1 (in-range (1+ x0) dim)) (x2 (in-range (1+ x1) dim))] (set! sol (+ (p5 x0) (p5 x1) (p5 x2)))
(when (vector-search sol all-y^5-x^5) ;; x0^5 + x1^5 + x2^5 = y^5 - x3^5 ??? (set! sol (append (list x0 x1 x2) (x-y sol))) ;; found (break #t))) ;; stop on first
→ (27 84 110 133 144) ;; time 2.8 sec
</lang>
ERRE
<lang ERRE>PROGRAM EULERO
CONST MAX=250
!$DOUBLE
FUNCTION POW5(X)
POW5=X*X*X*X*X
END FUNCTION
!$INCLUDE="PC.LIB"
BEGIN
CLS
FOR X0=1 TO MAX DO
FOR X1=1 TO X0 DO
FOR X2=1 TO X1 DO
FOR X3=1 TO X2 DO
LOCATE(3,1) PRINT(X0;X1;X2;X3)
SUM=POW5(X0)+POW5(X1)+POW5(X2)+POW5(X3)
S1=INT(SUM^0.2#+0.5#)
IF SUM=POW5(S1) THEN PRINT(X0,X1,X2,X3,S1) END IF
END FOR
END FOR
END FOR
END FOR
END PROGRAM</lang>
- Output:
133 110 84 27 144
F#
<lang fsharp> //Find 4 integers whose 5th powers sum to the fifth power of an integer (Quickly!) - Nigel Galloway: April 23rd., 2015 let G =
let GN = Array.init<float> 250 (fun n -> (float n)**5.0)
let rec gng (n, i, g, e) =
match (n, i, g, e) with
| (250,_,_,_) -> "No Solution Found"
| (_,250,_,_) -> gng (n+1, n+1, n+1, n+1)
| (_,_,250,_) -> gng (n, i+1, i+1, i+1)
| (_,_,_,250) -> gng (n, i, g+1, g+1)
| _ -> let l = GN.[n] + GN.[i] + GN.[g] + GN.[e]
match l with
| _ when l > GN.[249] -> gng(n,i,g+1,g+1)
| _ when l = round(l**0.2)**5.0 -> sprintf "%d**5 + %d**5 + %d**5 + %d**5 = %d**5" n i g e (int (l**0.2))
| _ -> gng(n,i,g,e+1)
gng (1, 1, 1, 1)
</lang>
- Output:
"27**5 + 84**5 + 110**5 + 133**5 = 144**5"
Fortran
<lang fortran>program sum_of_powers
implicit none
integer, parameter :: maxn = 249 integer, parameter :: dprec = selected_real_kind(15) integer :: i, x0, x1, x2, x3, y real(dprec) :: n(maxn), sumx
n = (/ (real(i, dprec)**5, i = 1, maxn) /)
outer: do x0 = 1, maxn
do x1 = 1, maxn
do x2 = 1, maxn
do x3 = 1, maxn
sumx = n(x0)+ n(x1)+ n(x2)+ n(x3)
y = 1
do while(y <= maxn .and. n(y) <= sumx)
if(n(y) == sumx) then
write(*,*) x0, x1, x2, x3, y
exit outer
end if
y = y + 1
end do
end do
end do
end do
end do outer
end program</lang>
- Output:
27 84 110 133 144
Go
<lang go>package main
import ( "fmt" "log" )
func main() { fmt.Println(eulerSum()) }
func eulerSum() (x0, x1, x2, x3, y int) { var pow5 [250]int for i := range pow5 { pow5[i] = i * i * i * i * i } for x0 = 4; x0 < len(pow5); x0++ { for x1 = 3; x1 < x0; x1++ { for x2 = 2; x2 < x1; x2++ { for x3 = 1; x3 < x2; x3++ { sum := pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3] for y = x0 + 1; y < len(pow5); y++ { if sum == pow5[y] { return } } } } } } log.Fatal("no solution") return }</lang>
- Output:
133 110 84 27 144
Haskell
<lang haskell>import Data.List import Data.List.Ordered
main :: IO () main = print $ head [(x0,x1,x2,x3,x4) |
-- choose x0, x1, x2, x3
-- so that 250 < x3 < x2 < x1 < x0
x3 <- [1..250-1],
x2 <- [1..x3-1],
x1 <- [1..x2-1],
x0 <- [1..x1-1],
let p5Sum = x0^5 + x1^5 + x2^5 + x3^5,
-- lazy evaluation of powers of 5
let p5List = [i^5|i <- [1..]],
-- is sum a power of 5 ?
member p5Sum p5List,
-- which power of 5 is sum ?
let Just x4 = elemIndex p5Sum p5List ]</lang>
- Output:
(27,84,110,133,144)
J
<lang J> (#~ (= <.)@((+/"1)&.:(^&5)))a.i.(#~ ] -:"1 /:~"1)@:([: ,/ ,"0 _1/)/ 4#,:}.250{.a. 27 84 110 133</lang>
Explanation:
There are 164059875 distinct possibilities to consider here (248 (+ %&(*/) ]) 1 2 3 4x) - so we need about 5GB memory to represent all the arguments using 64 bit integers. That's fairly trivial on some current laptops.
So, <lang J> a.i.(#~ ] -:"1 /:~"1)@:([: ,/ ,"0 _1/)/ 4#,:}.250{.a.</lang> finds all the possibilities for our four arguments and constrains them down to the distinct possibilities at each combining step. We use a character representation for our intermediate results, to limit memory overhead. Since this is a one-off problem, it's good enough to discard sequences which are not sorted.
Then, <lang J> (#~ (= <.)@((+/"1)&.:(^&5)))</lang> discards the cases we are not interested in. (It only keeps the case(s) where the fifth root of the sum of the fifth powers is an integer.)
Only one possibility remains.
JavaScript
<lang javascript>
var eulers_sum_of_powers = function (iMaxN) {
var aPow5 = []; var oPow5ToN = {};
for (var iP = 1; iP <= iMaxN; iP ++) { var iPow5 = Math.pow(iP, 5); aPow5.push(iPow5); oPow5ToN[iPow5] = iP; }
for (var i0 = 1; i0 <= iMaxN; i0 ++) { for (var i1 = 1; i1 <= i0; i1 ++) { for (var i2 = 1; i2 <= i1; i2 ++) { for (var i3 = 1; i3 <= i2; i3 ++) { var iPow5Sum = aPow5[i0] + aPow5[i1] + aPow5[i2] + aPow5[i3]; if (typeof oPow5ToN[iPow5Sum] != 'undefined') { return { i0: i0, i1: i1, i2: i2, i3: i3, iSum: oPow5ToN[iPow5Sum], }; } } } } }
};
var oResult = eulers_sum_of_powers(250);
console.log(oResult.i0 + '^5 + ' + oResult.i1 + '^5 + ' + oResult.i2 + '^5 + ' + oResult.i3 + '^5 = ' + oResult.iSum + '^5');
</lang>
Julia
<lang Julia> const lim = 250 const pwr = 5 const p = [i^pwr for i in 1:lim]
x = zeros(Int, pwr-1) y = 0
for a in combinations(1:lim, pwr-1)
b = searchsorted(p, sum(p[a])) 0 < length(b) || continue x = a y = b[1] break
end
if y == 0
println("No solution found for power = ", pwr, " and limit = ", lim, ".")
else
s = [@sprintf("%d^%d", i, pwr) for i in x]
s = join(s, " + ")
println("A solution is ", s, " = ", @sprintf("%d^%d", y, pwr), ".")
end </lang>
- Output:
A solution is 27^5 + 84^5 + 110^5 + 133^5 = 144^5.
Oforth
<lang Oforth>: eulerSum { | i j k l ip jp kp |
250 loop: i [
i 5 pow ->ip
i 1 + 250 for: j [
j 5 pow ip + ->jp
j 1 + 250 for: k [
k 5 pow jp + ->kp
k 1 + 250 for: l [
kp l 5 pow + 0.2 powf dup asInteger == ifTrue: [ [ i, j, k, l ] println ]
]
]
]
]
}</lang>
- Output:
>eulerSum [27, 84, 110, 133]
Pascal
slightly improved.Reducing calculation time by temporary sum and early break. <lang pascal>program Pot5Test; {$IFDEF FPC} {$MODE DELPHI}{$ELSE]{$APPTYPE CONSOLE}{$ENDIF} type
tTest = double;//UInt64;{ On linux 32Bit double is faster than Uint64 }
var
Pot5 : array[0..255] of tTest; res,tmpSum : tTest; x0,x1,x2,x3, y : NativeUint;//= Uint32 or 64 depending on OS xx-Bit i : byte;
BEGIN
For i := 1 to 255 do Pot5[i] := (i*i*i*i)*Uint64(i);
For x0 := 1 to 250-3 do
For x1 := x0+1 to 250-2 do
For x2 := x1+1 to 250-1 do
Begin
//set y here only, because pot5 is strong monoton growing,
//therefor the sum is strong monoton growing too.
y := x2+2;// aka x3+1
tmpSum := Pot5[x0]+Pot5[x1]+Pot5[x2];
For x3 := x2+1 to 250 do
Begin
res := tmpSum+Pot5[x3];
while (y< 250) AND (res > Pot5[y]) do
inc(y);
IF y > 250 then BREAK;
if res = Pot5[y] then
writeln(x0,'^5+',x1,'^5+',x2,'^5+',x3,'^5 = ',y,'^5');
end;
end;
END. </lang>
- output
27^5+84^5+110^5+133^5 = 144^5
real 0m1.091s {Uint64; Linux 32}real 0m0.761s {double; Linux 32}real 0m0.511s{Uint64; Linux 64}
Perl 6
<lang perl6>constant MAX = 250;
my %p5{Int}; my %sum2{Int};
for 1..MAX -> $i {
%p5{$i**5} = $i;
for 1..MAX -> $j {
%sum2{$i**5 + $j**5} = ($i, $j);
}
}
my @sk = %sum2.keys.sort; for %p5.keys.sort -> $p {
for @sk -> $s {
next if $p <= $s;
if %sum2{$p - $s} {
say (%sum2{$s}[],%sum2{$p-$s}[] X~ '⁵').join(' + ') ~ " = %p5{$p}⁵";
exit;
}
}
}</lang>
- Output:
84⁵ + 27⁵ + 133⁵ + 110⁵ = 144⁵
PHP
<lang php><?php
function eulers_sum_of_powers () { $max_n = 250; $pow_5 = array(); $pow_5_to_n = array(); for ($p = 1; $p <= $max_n; $p ++) { $pow5 = pow($p, 5); $pow_5 [$p] = $pow5; $pow_5_to_n[$pow5] = $p; } foreach ($pow_5 as $n_0 => $p_0) { foreach ($pow_5 as $n_1 => $p_1) { if ($n_0 < $n_1) continue; foreach ($pow_5 as $n_2 => $p_2) { if ($n_1 < $n_2) continue; foreach ($pow_5 as $n_3 => $p_3) { if ($n_2 < $n_3) continue; $pow_5_sum = $p_0 + $p_1 + $p_2 + $p_3; if (isset($pow_5_to_n[$pow_5_sum])) { return array($n_0, $n_1, $n_2, $n_3, $pow_5_to_n[$pow_5_sum]); } } } } } }
list($n_0, $n_1, $n_2, $n_3, $y) = eulers_sum_of_powers();
echo "$n_0^5 + $n_1^5 + $n_2^5 + $n_3^5 = $y^5";
?></lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
Python
<lang python>def eulers_sum_of_powers():
max_n = 250
pow_5 = [n**5 for n in range(max_n)]
pow5_to_n = {n**5: n for n in range(max_n)}
for x0 in range(1, max_n):
for x1 in range(1, x0):
for x2 in range(1, x1):
for x3 in range(1, x2):
pow_5_sum = sum(pow_5[i] for i in (x0, x1, x2, x3))
if pow_5_sum in pow5_to_n:
y = pow5_to_n[pow_5_sum]
return (x0, x1, x2, x3, y)
print("%i**5 + %i**5 + %i**5 + %i**5 == %i**5" % eulers_sum_of_powers())</lang>
- Output:
133**5 + 110**5 + 84**5 + 27**5 == 144**5
The above can be written as:
<lang python>from itertools import combinations
def eulers_sum_of_powers():
max_n = 250
pow_5 = [n**5 for n in range(max_n)]
pow5_to_n = {n**5: n for n in range(max_n)}
for x0, x1, x2, x3 in combinations(range(1, max_n), 4):
pow_5_sum = sum(pow_5[i] for i in (x0, x1, x2, x3))
if pow_5_sum in pow5_to_n:
y = pow5_to_n[pow_5_sum]
return (x0, x1, x2, x3, y)
print("%i**5 + %i**5 + %i**5 + %i**5 == %i**5" % eulers_sum_of_powers())</lang>
- Output:
27**5 + 84**5 + 110**5 + 133**5 == 144**5
It's much faster to cache and look up sums of two fifth powers, due to the small allowed range: <lang python>MAX = 250 p5, sum2 = {}, {}
for i in range(1, MAX): p5[i**5] = i for j in range(i, MAX): sum2[i**5 + j**5] = (i, j)
sk = sorted(sum2.keys()) for p in sorted(p5.keys()): for s in sk: if p <= s: break if p - s in sum2: print(p5[p], sum2[s] + sum2[p-s]) exit()</lang>
- Output:
144 (27, 84, 110, 133)
Racket
<lang scheme>#lang racket (define MAX 250) (define pow5 (make-vector MAX)) (for ([i (in-range 1 MAX)])
(vector-set! pow5 i (expt i 5)))
(define pow5s (list->set (vector->list pow5))) (let/ec break
(for* ([x0 (in-range 1 MAX)]
[x1 (in-range 1 x0)]
[x2 (in-range 1 x1)]
[x3 (in-range 1 x2)])
(define sum (+ (vector-ref pow5 x0)
(vector-ref pow5 x1)
(vector-ref pow5 x2)
(vector-ref pow5 x3)))
(when (set-member? pow5s sum)
(displayln (list x0 x1 x2 x3 (expt sum 1/5)))
(break))))</lang>
- Output:
(133 110 84 27 144.00000000000003)
REXX
Programming note: the 3rd argument can be specified which causes an attempt to find N solutions; the default is to find 1 solution.
The method used is:
- precompute all powers of five (within the confines of allowed integers)
- precompute all (positive) differences between two applicable 5thpowers
- see if any of the sums of any three 5th powers are equal to any of those (above) differences
- {thanks to the real nifty idea (↑↑↑) from userID G. Brougnard}
- see if the sum of any four 5th powers is equal to any 5th power
- {all of the above utilizes REXX's sparse stemmed array hashing which eliminates the need for sorting.}
By implementing (userID) G. Brougnard's idea of differences of two 5th powers, the time used for computation was cut by over a factor of seventy.
In essence, the formula being solved is: aⁿ + bⁿ + cⁿ == xⁿ ─ dⁿ
which lends itself to algorithm optimization by (only) having to pre-compute all possible differences between two integer powers of five, and then (only) summing three integer powers of five.
<lang rexx>/*REXX program finds unique positive integer solution(s) for the equation: ─────
────────────────────────────────────── aⁿ + bⁿ + cⁿ + dⁿ == xⁿ where n=5*/
parse arg L H N . /*get optional LOW, HIGH, #solutions.*/
if L== | L==',' then L= 0 + 1 /*Not specified? Then use the default.*/
if H== | H==',' then H=250 - 1 /* " " " " " " */
if N== | N==',' then N= 1 /* " " " " " " */
numeric digits 1000 /*be able to handle the next expression*/
numeric digits max(9,length(4*(H**5))) /* " " " " integer 5th powers.*/
h1=H-1; h2=H-2; h3=H-3 /*calculate the upper DO loop limits.*/
say center(' 'strip(translate(sourceLine(2), ,'* /'))' ', 79, ' ') /*title.*/
/* [↓] define values of 5th powers. */
!.=0; do pow=L to H; @.pow=pow**5; _=@.pow; !._=1; #._=pow; end ?.=0; do j=L+3 to H-1; do k=j+1 to H; _=@.k-@.j; ?._=1; end; end
/* [↓] define 5th power differences. */
- =0 /*#: is the number of solutions found.*/
do a=L to h3; s0= @.a /*traipse through possible A values. */
do b=a+1 to h2; s1=s0+@.b /* " " " B " */
do c=b+1 to h1; s2=s1+@.c /* " " " C " */
if ?.s2 then do d=c+1 to H; s3=s2+@.d /*find appropriate solution.*/
if !.s3 then call results /*A solution? Then show it.*/
end /*d*/ /* [↑] !.S3 is a boolean.*/
end /*c*/
end /*b*/
end /*a*/
if #==0 then say "Didn't find a solution."; signal done /*────────────────────────────────────────────────────────────────────────────*/ results: _=left(,5); #=#+1 /*_ is spacing between values; bump #.*/ say _ 'solution' #":" _ 'a='a _ "b="||b _ 'c='c _ "d="d _ 'x='#.s3 if #<N then return /*return, keep searching for more sols.*/ done: exit # /*stick a fork in it, we're all done. */</lang> output when using the default inputs:
─────────────────── aⁿ + bⁿ + cⁿ + dⁿ == xⁿ where n=5 ───────────────────
solution 1: a=27 b=84 c=110 d=133 x=144
Ruby
Brute force: <lang ruby>power5 = (1..250).each_with_object({}){|i,h| h[i**5]=i} result = power5.keys.repeated_combination(4).select{|a| power5[a.inject(:+)]} puts result.map{|a| a.map{|i| "#{power5[i]}**5"}.join(' + ') + " = #{power5[a.inject(:+)]}**5"}</lang>
- Output:
27**5 + 84**5 + 110**5 + 133**5 = 144**5
Faster version:
<lang ruby>p5, sum2, max = {}, {}, 250 (1..max).each do |i|
p5[i**5] = i
(i..max).each{|j| sum2[i**5 + j**5] = [i,j]}
end
result = {} sk = sum2.keys.sort p5.keys.sort.each do |p|
sk.each do |s| break if p <= s result[(sum2[s] + sum2[p-s]).sort] = p5[p] if sum2[p - s] end
end result.each{|k,v| puts k.map{|i| "#{i}**5"}.join(' + ') + " = #{v}**5"}</lang> The output is the same above.
Rust
<lang rust>const MAX_N : u64 = 250;
fn eulers_sum_of_powers() -> (usize, usize, usize, usize, usize) {
let pow5: Vec<u64> = (0..MAX_N).map(|i| i.pow(5)).collect(); let pow5_to_n = |pow| pow5.binary_search(&pow);
for x0 in 1..MAX_N as usize {
for x1 in 1..x0 {
for x2 in 1..x1 {
for x3 in 1..x2 {
let pow_sum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3];
if let Ok(n) = pow5_to_n(pow_sum) {
return (x0, x1, x2, x3, n)
}
}
}
}
}
panic!();
}
fn main() { let (x0, x1, x2, x3, y) = eulers_sum_of_powers(); println!("{}^5 + {}^5 + {}^5 + {}^5 == {}^5", x0, x1, x2, x3, y) }</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 == 144^5
VBScript
<lang vb>Max=250
For X0=1 To Max For X1=1 To X0 For X2=1 To X1 For X3=1 To X2 Sum=fnP5(X0)+fnP5(X1)+fnP5(X2)+fnP5(X3) S1=Int(Sum^0.2) If Sum=fnP5(S1) Then WScript.StdOut.Write X0 & " " & X1 & " " & X2 & " " & X3 & " " & S1 WScript.Quit End If Next Next Next Next
Function fnP5(n) fnP5 = n ^ 5 End Function</lang>
- Output:
133 110 84 27 144
zkl
Uses two look up tables for efficiency. Counts from 0 for ease of coding. <lang zkl>pow5s:=[1..249].apply("pow",5); // (1^5, 2^5, 3^5 .. 249^5) pow5r:=pow5s.enumerate().apply("reverse").toDictionary(); // [144^5:144, ...] foreach x0,x1,x2,x3 in (249,x0,x1,x2){
sum:=pow5s[x0] + pow5s[x1] + pow5s[x2] + pow5s[x3];
if(pow5r.holds(sum))
println("%d^5 + %d^5 + %d^5 + %d^5 = %d^5"
.fmt(x3+1,x2+1,x1+1,x0+1,pow5r[sum]+1));
}</lang>
- Output:
27^5 + 84^5 + 110^5 + 133^5 = 144^5
Using the Python technique of caching double sums results in a 5x speed up [to the first/only solution]; actually the speed up is close to 25x but creating the caches dominates the runtime to the first solution.
<lang zkl>p5,sum2:=D(),D(); foreach i in ([1..249]){
p5[i.pow(5)]=i;
foreach j in ([i..249]){ sum2[i.pow(5) + j.pow(5)]=T(i,j) } // 31,125 keys
}
sk:=sum2.keys.apply("toInt").copy().sort(); // RW list sorts faster than a RO one foreach p,s in (p5.keys.apply("toInt"),sk){
if(p<=s) break;
if(sum2.holds(p - s)){
println("%d^5 + %d^5 + %d^5 + %d^5 = %d^5"
.fmt(sum2[s].xplode(),sum2[p - s].xplode(),p5[p]));
break(2); // or get permutations
}
}</lang> Note: dictionary keys are always strings and copying a read only list creates a read write list.
- Output:
27^5 + 84^5 + 110^5 + 133^5 = 144^5