# Euler's sum of powers conjecture

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.

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).

## 360 Assembly

This program could have been run in 1964. Here, for maximum compatibility, we use only the basic 360 instruction set. Macro instruction XPRNT can be replaced by a WTO. <lang 360asm> EULERCO CSECT

```        USING  EULERCO,R13
B      80(R15)
DC     17F'0'
DC     CL8'EULERCO'
STM    R14,R12,12(R13)
ST     R13,4(R15)
ST     R15,8(R13)
LR     R13,R15
ZAP    X1,=P'1'
```

LOOPX1 ZAP PT,MAXN do x1=1 to maxn-4

```        SP     PT,=P'4'
CP     X1,PT
BH     ELOOPX1
ZAP    PT,X1
AP     PT,=P'1'
ZAP    X2,PT
```

LOOPX2 ZAP PT,MAXN do x2=x1+1 to maxn-3

```        SP     PT,=P'3'
CP     X2,PT
BH     ELOOPX2
ZAP    PT,X2
AP     PT,=P'1'
ZAP    X3,PT
```

LOOPX3 ZAP PT,MAXN do x3=x2+1 to maxn-2

```        SP     PT,=P'2'
CP     X3,PT
BH     ELOOPX3
ZAP    PT,X3
AP     PT,=P'1'
ZAP    X4,PT
```

LOOPX4 ZAP PT,MAXN do x4=x3+1 to maxn-1

```        SP     PT,=P'1'
CP     X4,PT
BH     ELOOPX4
ZAP    PT,X4
AP     PT,=P'1'
ZAP    X5,PT              x5=x4+1
ZAP    SUMX,=P'0'         sumx=0
ZAP    PT,X1              x1
BAL    R14,POWER5
AP     SUMX,PT
ZAP    PT,X2              x2
BAL    R14,POWER5
AP     SUMX,PT
ZAP    PT,X3              x3
BAL    R14,POWER5
AP     SUMX,PT
ZAP    PT,X4              x4
BAL    R14,POWER5
AP     SUMX,PT            sumx=x1**5+x2**5+x3**5+x4**5
ZAP    PT,X5              x5
BAL    R14,POWER5
ZAP    VALX,PT            valx=x5**5
```

LOOPX5 CP X5,MAXN while x5<=maxn & valx<=sumx

```        BH     ELOOPX5
CP     VALX,SUMX
BH     ELOOPX5
CP     VALX,SUMX          if valx=sumx
BNE    NOTEQUAL
MVI    BUF,C' '
MVC    BUF+1(79),BUF      clear buffer
ED     WC,X1              x1
MVC    BUF+0(8),WC+8
ED     WC,X2              x2
MVC    BUF+8(8),WC+8
ED     WC,X3              x3
MVC    BUF+16(8),WC+8
ED     WC,X4              x4
MVC    BUF+24(8),WC+8
ED     WC,X5              x5
MVC    BUF+32(8),WC+8
XPRNT  BUF,80             output x1,x2,x3,x4,x5
B      ELOOPX1
```

NOTEQUAL ZAP PT,X5

```        AP     PT,=P'1'
ZAP    X5,PT              x5=x5+1
ZAP    PT,X5
BAL    R14,POWER5
ZAP    VALX,PT            valx=x5**5
B      LOOPX5
```

ELOOPX5 AP X4,=P'1'

```        B      LOOPX4
```

ELOOPX4 AP X3,=P'1'

```        B      LOOPX3
```

ELOOPX3 AP X2,=P'1'

```        B      LOOPX2
```

ELOOPX2 AP X1,=P'1'

```        B      LOOPX1
```

ELOOPX1 L R13,4(0,R13)

```        LM     R14,R12,12(R13)
XR     R15,R15
BR     R14
```

POWER5 ZAP PQ,PT ^1

```        MP     PQ,PT              ^2
MP     PQ,PT              ^3
MP     PQ,PT              ^4
MP     PQ,PT              ^5
ZAP    PT,PQ
BR     R14
```

MAXN DC PL8'250' X1 DS PL8 X2 DS PL8 X3 DS PL8 X4 DS PL8 X5 DS PL8 SUMX DS PL8 VALX DS PL8 PT DS PL8 PQ DS PL8 WC DS CL17 MASK DC X'40',13X'20',X'212060' CL17 BUF DS CL80

```        YREGS
END
```

</lang>

Output:
```      27      84     110     133     144
```

procedure Sum_Of_Powers is

```  type Base is range 0 .. 250; -- A, B, C, D and Y are in that range
type Num is range 0 .. 4*(250**5); -- (A**5 + ... + D**5) is in that range
subtype Fit is Num range 0 .. 250**5; -- Y**5 is in that range

Modulus: constant Num := 254;
type Modular is mod Modulus;

type Result_Type is array(1..5) of Base; -- this will hold A,B,C,D and Y

type Y_Type is array(Modular) of Base;
type Y_Sum_Type is array(Modular) of Fit;
```
```  Y_Sum: Y_Sum_Type := (others => 0);
Y: Y_Type := (others => 0);
-- for I in 0 .. 250, we set Y_Sum(I**5 mod Modulus) := I**5
--                       and Y(I**5 mod Modulus) := I
-- Modulus has been chosen to avoid collisions on (I**5 mod Modulus)
-- later we will compute Sum_ABCD := A**5 + B**5 + C**5 + D**5
-- and check if Y_Sum(Sum_ABCD mod modulus) = Sum_ABCD

function Compute_Coefficients return Result_Type is

Sum_A: Fit;
Sum_AB, Sum_ABC, Sum_ABCD: Num;
Short: Modular;

begin
for A in Base(0) .. 246 loop
Sum_A := Num(A) ** 5;
for B in A .. 247 loop
Sum_AB := Sum_A + (Num(B) ** 5);
for C in Base'Max(B,1) .. 248 loop -- if A=B=0 then skip C=0
Sum_ABC := Sum_AB + (Num(C) ** 5);
for D in C .. 249 loop
Sum_ABCD := Sum_ABC + (Num(D) ** 5);
Short    := Modular(Sum_ABCD mod Modulus);
if Y_Sum(Short) = Sum_ABCD then
return A & B & C & D & Y(Short);
end if;
end loop;
end loop;
end loop;
end loop;
return 0 & 0 & 0 & 0 & 0;
end Compute_Coefficients;
```
```  Tmp: Fit;
ABCD_Y: Result_Type;
```

begin -- main program

```  -- initialize Y_Sum and Y
for I in Base(0) .. 250 loop
Tmp := Num(I)**5;
if Y_Sum(Modular(Tmp mod Modulus)) /= 0 then
raise Program_Error with "Collision: Change Modulus and recompile!";
else
Y_Sum(Modular(Tmp mod Modulus)) := Tmp;
Y(Modular(Tmp mod Modulus)) := I;
end if;
end loop;

-- search for a solution (A, B, C, D, Y)
ABCD_Y := Compute_Coefficients;
```
```  -- output result
for Number of ABCD_Y loop
end loop;

```

end Sum_Of_Powers;</lang>

Output:
``` 27 84 110 133 144
```

## ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.win32

<lang algol68># max number will be the highest integer we will consider # INT max number = 250;

1. Construct a table of the fifth powers of 1 : max number #

[ max number ]LONG INT fifth; FOR i TO max number DO

```   LONG INT i2 =  i * i;
fifth[ i ] := i2 * i2 * i
```

OD;

1. find the first a, b, c, d, e such that a^5 + b^5 + c^5 + d^5 = e^5 #
2. as the fifth powers are in order, we can use a binary search to determine #
3. whether the value is in the table #

BOOL found := FALSE; FOR a TO max number WHILE NOT found DO

```   FOR b FROM a TO max number WHILE NOT found DO
LONG INT sum   = fifth[a] + fifth[b] + fifth[c] + fifth[d];
INT      low  := d;
INT      high := max number;
WHILE low < high
DO
INT e := ( low + high ) OVER 2;
IF fifth[ e ] = sum
THEN
# the value at e is a fifth power                    #
found := TRUE;
print( ( ( whole( a, 0 ) + "^5 + " + whole( b, 0 ) + "^5 + "
+ whole( c, 0 ) + "^5 + " + whole( d, 0 ) + "^5 = "
+ whole( e, 0 ) + "^5"
)
, newline
)
)
ELIF sum < fifth[ e ]
THEN high := e - 1
ELSE low  := e + 1
FI
OD
OD
OD
OD
```

OD</lang> Output:

```27^5 + 84^5 + 110^5 + 133^5 = 144^5
```

## 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=e0+30k, for some starting value e0, and any k. Also, we follow the constraints 1<=a<b<c<d<e<N in the main loop. <lang c>// Alexander Maximov, July 2nd, 2015

1. include <stdio.h>
2. include <time.h>

typedef long long mylong;

void compute(int N, char find_only_one_solution) { const int M = 30; /* x^5 == x modulo M=2*3*5 */ int a, b, c, d, e; mylong s, t, max, *p5 = (mylong*)malloc(sizeof(mylong)*(N+M));

for(s=0; s < N; ++s) p5[s] = s * s, p5[s] *= p5[s] * s; for(max = p5[N - 1]; s < (N + M); p5[s++] = max + 1);

for(a = 1; a < N; ++a) for(b = a + 1; b < N; ++b) for(c = b + 1; c < N; ++c) for(d = c + 1, e = d + ((t = p5[a] + p5[b] + p5[c]) % M); ((s = t + p5[d]) <= max); ++d, ++e) { for(e -= M; p5[e + M] <= 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); if(find_only_one_solution) goto onexit; } } onexit: free(p5); }

int main(void) { int tm = clock(); compute(250, 0); printf("time=%d milliseconds\r\n", (int)((clock() - tm) * 1000 / CLOCKS_PER_SEC)); return 0; }</lang>

Output:

The fair way to measure the speed of the code above is to measure it's run time to find all possible solutions to the problem, given N (and not just a single solution, since then the time may depend on the way and the order we organize for-loops).

```27 84 110 133 144
time=235 milliseconds
```

Another test with N=1000 produces the following results:

```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
time=65743 milliseconds
```

PS: The solution for C++ provided below is actually quite good in its design idea behind. However, with all proposed tricks to speed up, the measurements for C++ solution for N=1000 showed the execution time 81447ms (+23%) on the same environment as above for C solution (same machine, same compiler, 64-bit platform). The reason that C++ solution is a bit slower is, perhaps, the fact that the inner loops over rs have complexity ~N/2 steps in average, while with the modulo 30 trick that complexity can be reduced down to ~N/60 steps, although one "expensive" extra %-operation is still needed.

## C++

The simplest brute-force find is already reasonably quick: <lang cpp>#include <iostream>

1. include <cmath>
2. include <set>

using namespace std;

bool find() { const auto MAX = 250; vector<double> pow5(MAX); for (auto i = 1; i < MAX; i++) pow5[i] = (double)i * i * i * i * 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 (binary_search(pow5.begin(), pow5.end(), sum)) { cout << x0 << " " << x1 << " " << x2 << " " << x3 << " " << pow(sum, 1.0 / 5.0) << endl; return true; } } } } } // not found return false; }

int main(void) { int tm = clock(); if (!find()) cout << "Nothing found!\n"; printf("time=%d milliseconds\r\n", (int)((clock() - tm) * 1000 / CLOCKS_PER_SEC)); return 0; } </lang>

Output:
```133 110 84 27 144
time=234 milliseconds
```

We can accelerate this further by creating a parallel std::set<double> p5s containing the elements of the std::vector pow5, and using it to replace the call to std::binary_search: <lang> set<double> pow5s; for (auto i = 1; i < MAX; i++) { pow5[i] = (double)i * i * i * i * i; pow5s.insert(pow5[i]); } //...

```       if (pow5s.find(sum) != pow5s.end())
```

</lang> This reduces the timing to 125 ms on the same hardware.

A different, more effective optimization is to note that each successive sum is close to the previous one, and use a bidirectional linear search with memory. <lang> bool find() { const auto MAX = 250; vector<double> pow5(MAX); for (auto i = 1; i < MAX; i++) pow5[i] = (double)i * i * i * i * i; auto rs = 5; 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]; while (rs > 0 && pow5[rs] > sum) --rs; while (rs < MAX - 1 && pow5[rs] < sum) ++rs; if (pow5[rs] == sum) { cout << x0 << " " << x1 << " " << x2 << " " << x3 << " " << pow(sum, 1.0 / 5.0) << endl; return true; } } } } } // not found return false; } </lang> This reduces the timing to 47 ms. We could equally well replace rs with an iterator inside pow5; the timing is unaffected.

Finally, we can set rs to a lower value in the loop over x2, and then never have to decrement it during the loop over x3: <lang> for (auto x2 = 1; x2 < x1; x2++) { auto s2 = pow5[x0] + pow5[x1] + pow5[x2]; while (rs > 0 && pow5[rs] > s2) --rs; for (auto x3 = 1; x3 < x2; x3++) { auto sum = s2 + pow5[x3]; while (rs < MAX - 1 && pow5[rs] < sum) ++rs; if (pow5[rs] == sum) </lang> This approximately halves the execution time, which is now too short to reliably measure.

## 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

Works with: Fortran version 95 and later

<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

Translation of: Python

<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
```

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.

## Java

Translation of: ALGOL 68

Tested with Java 6. <lang java>public class eulerSopConjecture {

```   static final int    MAX_NUMBER = 250;
```
```   public static void main( String[] args )
{
boolean found = false;
long[]  fifth = new long[ MAX_NUMBER ];
```
```       for( int i = 1; i <= MAX_NUMBER; i ++ )
{
long i2 =  i * i;
fifth[ i - 1 ] = i2 * i2 * i;
} // for i
```
```       for( int a = 0; a < MAX_NUMBER && ! found ; a ++ )
{
for( int b = a; b < MAX_NUMBER && ! found ; b ++ )
{
for( int c = b; c < MAX_NUMBER && ! found ; c ++ )
{
for( int d = c; d < MAX_NUMBER && ! found ; d ++ )
{
long sum  = fifth[a] + fifth[b] + fifth[c] + fifth[d];
int  e = java.util.Arrays.binarySearch( fifth, sum );
found  = ( e >= 0 );
if( found )
{
// the value at e is a fifth power
System.out.print( (a+1) + "^5 + "
+ (b+1) + "^5 + "
+ (c+1) + "^5 + "
+ (d+1) + "^5 = "
+ (e+1) + "^5"
);
} // if found;;
} // for d
} // for c
} // for b
} // for a
} // main
```

} // eulerSopConjecture</lang> Output:

```27^5 + 84^5 + 110^5 + 133^5 = 144^5
```

## JavaScript

 This example does not show the output mentioned in the task description on this page (or a page linked to from here). Please ensure that it meets all task requirements and remove this message. Note that phrases in task descriptions such as "print and display" and "print and show" for example, indicate that (reasonable length) output be a part of a language's solution.

<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>

## jq

Works with: jq version 1.4

This version finds all non-decreasing solutions within the specified bounds, using a brute-force but not entirely blind approach. <lang jq># Search for y in 1 .. maxn (inclusive) for a solution to SIGMA (xi ^ 5) = y^5

1. and for each solution with x0<=x1<=...<x3, print [x0, x1, x3, x3, y]

def sum_of_powers_conjecture(maxn):

``` def p5: . as \$in | (.*.) | ((.*.) * \$in);
def fifth: log / 5 | exp;
```
``` # return the fifth root if . is a power of 5
def integral_fifth_root: fifth | if . == floor then . else false end;
```
``` (maxn | p5) as \$uber
| range(1; maxn) as \$x0
| (\$x0 | p5) as \$s0
| if \$s0 < \$uber then range(\$x0; (\$uber - \$s0 | fifth) + 1) as \$x1
| (\$s0 + (\$x1 | p5)) as \$s1
| if \$s1 < \$uber then range(\$x1; (\$uber - \$s1 | fifth) + 1) as \$x2
| (\$s1 + (\$x2 | p5)) as \$s2
| if \$s2 < \$uber then range(\$x2; (\$uber - \$s2 | fifth) + 1) as \$x3
| (\$s2 + (\$x3 | p5)) as \$sumx
```

| (\$sumx | integral_fifth_root) | if . then [\$x0,\$x1,\$x2,\$x3,.] else empty end else empty end

```     else empty
end
else empty
end ;</lang>
```

Output:

<lang sh>\$ jq -c -n -f Euler_sum_of_powers_conjecture_fifth_root.jq [27,84,110,133,144]</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

Works with: Free 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

Translation of: Python

<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

Translation of: Python

<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:

Works with: Python version 2.6+

<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

 This example is incomplete. Integer solutions required. Please ensure that it meets all task requirements and remove this message.
Translation of: C++

<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 new 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 any two applicable integer powers of five (which happens to be 30,135 unique differences), and then (only) summing any applicable 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.  */
```
1. =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:

Translation of: Python

<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

Translation of: ERRE

<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.

Translation of: Python

<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`