Euler's sum of powers conjecture: Difference between revisions
→{{header|ZX Spectrum Basic}}: no need to follow calculations for the QL |
|||
Line 4,759:
5 DIM k(29): DIM q(249)
10 FOR i=4 TO 249: LET q(i)=LN i : NEXT i
11 REM enhancements for the much expanded Spectrum Next: DIM p(
12 REM FOR j=
14 PRINT "slide rule ready"
15 FOR i=0 TO 9: LET k(i)=240+ i : NEXT i
Line 4,779:
45 LET s=EXP((q(y)-q(m))*5)+ s
47 LET s=EXP((q(z)-q(m))*5)+ s
50 IF s<>1 THEN GO TO 80
52 LET a=FN f(w*w,m) : LET a=FN f(a*
53 LET b=FN f(x*x,m) : LET b=FN f(b*
55 LET c=FN f(y*y,m) : LET c=FN f(c*
57 LET d=FN f(z*z,m) : LET d=FN f(d*
60 LET u=FN f((a+b+c+d),m)
65 IF u THEN GO TO 80
|
Revision as of 00:28, 25 December 2021
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 two hundred 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
In 1966, Leon J. Lander and Thomas R. Parkin 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 for:
x05 + x15 + x25 + x35 == y5
Where all xi
's and y
are distinct integers between 0 and 250 (exclusive).
Show an answer here.
- Related tasks
11l
<lang 11l>F eulers_sum_of_powers()
V max_n = 250 V pow_5 = (0 .< max_n).map(n -> Int64(n) ^ 5) V pow5_to_n = Dict(0 .< max_n, n -> (Int64(n) ^ 5, n))
L(x0) 1 .< max_n L(x1) 1 .< x0 L(x2) 1 .< x1 L(x3) 1 .< x2 V pow_5_sum = pow_5[x0] + pow_5[x1] + pow_5[x2] + pow_5[x3] I pow_5_sum C pow5_to_n V y = pow5_to_n[pow_5_sum] R (x0, x1, x2, x3, y)
V r = eulers_sum_of_powers() print(‘#.^5 + #.^5 + #.^5 + #.^5 = #.^5’.format(r[0], r[1], r[2], r[3], r[4]))</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
360 Assembly
In the program we do not user System/360 integers (31 bits) unable to handle the problem, but System/360 packed decimal (15 digits). 250^5 needs 12 digits.
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 MVC WC,MASK ED WC,X1 x1 MVC BUF+0(8),WC+8 MVC WC,MASK ED WC,X2 x2 MVC BUF+8(8),WC+8 MVC WC,MASK ED WC,X3 x3 MVC BUF+16(8),WC+8 MVC WC,MASK ED WC,X4 x4 MVC BUF+24(8),WC+8 MVC WC,MASK 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'249' 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
6502 Assembly
This is long enough as it is, so the code here is just the main task body. Details like the BASIC loader (for a C-64, which is what I ran this on) and the output routines (including the very tedious conversion of a 40-bit integer to decimal) were moved into external include files. Also in an external file is the prebuilt table of the first 250 fifth powers ... actually, just for ease of referencing, it's the first 256, 0...255, in five tables each holding one byte of the value.
<lang assembly>; Prove Euler's sum of powers conjecture false by finding
- positive a,b,c,d,e such that a⁵+b⁵+c⁵+d⁵=e⁵.
- we're only looking for the first counterexample, which occurs with all
- integers less than this value
max_value = $fa ; decimal 250
- this header turns our code into a LOADable and RUNnable BASIC program
.include "basic_header.s"
- this contains the first 256 integers to the power of 5 broken up into
- 5 tables of one-byte values (power5byte0 with the LSBs through
- power5byte4 with the MSBs)
.include "power5table.s"
- this defines subroutines and macros for printing messages to
- the console, including `puts` for printing out a NUL-terminated string,
- puthex to display a one-byte value in hexadecimal, and putdec through
- putdec5 to display an N-byte value in decimal
.include "output.s"
- label strings for the result output
.feature string_escapes success: .asciiz "\r\rFOUND EXAMPLE:\r\r" between: .asciiz "^5 + " penult: .asciiz "^5 = " eqend: .asciiz "^5\r\r(SUM IS "
- the BASIC loader program prints the elapsed time at the end, so we include a
- label for that, too
tilabel: .asciiz ")\r\rTIME:"
- ZP locations to store the integers to try
x0 = $f7 x1 = x0 + 1 x2 = x0 + 2 x3 = x0 + 3
- we use binary search to find integer roots; current bounds go here
low = x0 + 4 hi = x0 + 5
- sum of powers of current candidate integers
sum: .res 5
- when we find a sum with an integer 5th root, we put it here
x4: .res 1
main: ; loop for x0 from 1 to max_value
ldx #01 stx x0
loop0: ; loop for x1 from x0+1 to max_value
ldx x0 inx stx x1
loop1: ; loop for x2 from x1+1 to max_value
ldx x1 inx stx x2
loop2: ; loop for x3 from x2+1 to max_value
ldx x2 inx stx x3
loop3: ; add up the fifth powers of the four numbers
; initialize to 0 lda #00 sta sum sta sum+1 sta sum+2 sta sum+3 sta sum+4
; we use indexed addressing, taking advantage of the fact that the xn's ; are consecutive, so x0,1 = x1, etc. ldy #0
addloop:
ldx x0,y lda sum clc adc power5byte0,x sta sum lda sum+1 adc power5byte1,x sta sum+1 lda sum+2 adc power5byte2,x sta sum+2 lda sum+3 adc power5byte3,x sta sum+3 lda sum+4 adc power5byte4,x sta sum+4 iny cpy #4 bcc addloop ; now sum := x₀⁵+x₁⁵+x₂⁵+x₃⁵ ; set initial bounds for binary search ldx x3 inx stx low ldx #max_value dex stx hi
binsearch:
; compute midpoint lda low cmp hi beq notdone bcs done_search
notdone:
ldx #0 clc adc hi
; now a + carry bit = low+hi; rotating right will get the midpoint ror ; compare square of midpoint to sum tax lda sum+4 cmp power5byte4,x bne notyet lda sum+3 cmp power5byte3,x bne notyet lda sum+2 cmp power5byte2,x bne notyet lda sum+1 cmp power5byte1,x bne notyet lda sum cmp power5byte0,x beq found
notyet:
bcc sum_lt_guess inx stx low bne endbin beq endbin
sum_lt_guess:
dex stx hi
endbin:
bne binsearch beq binsearch
done_search:
inc x3 lda x3 cmp #max_value bcs end_loop3 jmp loop3
end_loop3:
inc x2 lda x2 cmp #max_value bcs end_loop2 jmp loop2
end_loop2:
inc x1 lda x1 cmp #max_value bcs end_loop1 jmp loop1
end_loop1:
inc x0 lda x0 cmp #max_value bcs end_loop0 jmp loop0
end_loop0:
; should never get here, means we didn't find an example. brk
found: stx x4
puts success putdec x0 ldy #1
ploop: puts between
putdec {x0,y} iny cpy #4 bcc ploop puts penult putdec {x0,y} puts eqend putdec5 sum puts tilabel rts</lang>
- Output:
**** COMMODORE 64 BASIC V2 **** 64K RAM SYSTEM 38911 BASIC BYTES FREE READY. LOAD"ESOP",8,1 SEARCHING FOR ESOP LOADING READY. RUN: FOUND EXAMPLE: 27^5 + 84^5 + 110^5 + 133^5 = 144^5 (SUM IS 61917364224) TIME:095617 READY.
... almost ten hours, but it did find it!
Ada
<lang Ada>with Ada.Text_IO;
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 Ada.Text_IO.Put(Base'Image(Number)); end loop; Ada.Text_IO.New_Line;
end Sum_Of_Powers;</lang>
- Output:
27 84 110 133 144
ALGOL 68
<lang algol68># max number will be the highest integer we will consider # INT max number = 250;
- 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;
- find the first a, b, c, d, e such that a^5 + b^5 + c^5 + d^5 = e^5 #
- as the fifth powers are in order, we can use a binary search to determine #
- 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 FOR c FROM b TO max number WHILE NOT found DO FOR d FROM c 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 AND NOT found 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
ALGOL W
As suggested by the REXX solution, we find a solution to a^5 + b^5 + c^5 = e^5 - d^5 which results in a significant reduction in run time.
Algol W integers are 32-bit only, so we simulate the necessary 12 digit arithmetic with pairs of integers.
<lang algolw>begin
% find a, b, c, d, e such that a^5 + b^5 + c^5 + d^5 = e^5 % % where 1 <= a <= b <= c <= d <= e <= 250 % % we solve this using the equivalent equation a^5 + b^5 + c^5 = e^5 - d^5 % % 250^5 is 976 562 500 000 - too large for a 32 bit number so we will use pairs of % % integers and constrain their values to be in 0..1 000 000 % % Note only positive numbers are needed % integer MAX_NUMBER, MAX_V; MAX_NUMBER := 250; MAX_V := 1000000; begin % quick sorts the fifth power differences table % procedure quickSort5 ( integer value lb, ub ) ; if ub > lb then begin % more than one element, so must sort % integer left, right, pivot, pivotLo, pivotHi; left := lb; right := ub; % choosing the middle element of the array as the pivot % pivot := left + ( ( ( right + 1 ) - left ) div 2 ); pivotLo := loD( pivot ); pivotHi := hiD( pivot ); while begin while left <= ub and begin integer cmp; cmp := hiD( left ) - pivotHi; if cmp = 0 then cmp := loD( left ) - pivotLo; cmp < 0 end do left := left + 1; while right >= lb and begin integer cmp; cmp := hiD( right ) - pivotHi; if cmp = 0 then cmp := loD( right ) - pivotLo; cmp > 0 end do right := right - 1; left <= right end do begin integer swapLo, swapHi, swapD, swapE; swapLo := loD( left ); swapHi := hiD( left ); swapD := Dd( left ); swapE := De( left ); loD( left ) := loD( right ); hiD( left ) := hiD( right ); Dd( left ) := Dd( right ); De( left ) := De( right ); loD( right ) := swapLo; hiD( right ) := swapHi; Dd( right ) := swapD; De( right ) := swapE; left := left + 1; right := right - 1 end while_left_le_right ; quickSort5( lb, right ); quickSort5( left, ub ) end quickSort5 ; % table of fifth powers % integer array lo5, hi5 ( 1 :: MAX_NUMBER ); % table if differences between fifth powers % integer array loD, hiD, De, Dd ( 1 :: MAX_NUMBER * MAX_NUMBER ); integer dUsed, dPos; % compute fifth powers % for i := 1 until MAX_NUMBER do begin lo5( i ) := i * i; hi5( i ) := 0; for p := 3 until 5 do begin integer carry; lo5( i ) := lo5( i ) * i; carry := lo5( i ) div MAX_V; lo5( i ) := lo5( i ) rem MAX_V; hi5( i ) := hi5( i ) * i; hi5( i ) := hi5( i ) + carry end for_p end for_i ; % compute the differences between fifth powers e^5 - d^5, 1 <= d < e <= MAX_NUMBER % dUsed := 0; for e := 2 until MAX_NUMBER do begin for d := 1 until e - 1 do begin dUsed := dUsed + 1; De( dUsed ) := e; Dd( dUsed ) := d; loD( dUsed ) := lo5( e ) - lo5( d ); hiD( dUsed ) := hi5( e ) - hi5( d ); if loD( dUsed ) < 0 then begin loD( dUsed ) := loD( dUsed ) + MAX_V; hiD( dUsed ) := hiD( dUsed ) - 1 end if_need_to_borrow end for_d end for_e; % sort the fifth power differences % quickSort5( 1, dUsed ); % attempt to find a^5 + b^5 + c^5 = e^5 - d^5 % for a := 1 until MAX_NUMBER do begin integer loA, hiA; loA := lo5( a ); hiA := hi5( a ); for b := a until MAX_NUMBER do begin integer loB, hiB; loB := lo5( b ); hiB := hi5( b ); for c := b until MAX_NUMBER do begin integer low, high, loSum, hiSum; loSum := loA + loB + lo5( c ); hiSum := ( loSum div MAX_V ) + hiA + hiB + hi5( c ); loSum := loSum rem MAX_V; % look for hiSum,loSum in hiD,loD % low := 1; high := dUsed; while low < high do begin integer mid, cmp; mid := ( low + high ) div 2; cmp := hiD( mid ) - hiSum; if cmp = 0 then cmp := loD( mid ) - loSum; if cmp = 0 then begin % the value at mid is the difference of two fifth powers % write( i_w := 1, s_w := 0 , a, "^5 + ", b, "^5 + ", c, "^5 + " , Dd( mid ), "^5 = ", De( mid ), "^5" ); go to found end else if cmp > 0 then high := mid - 1 else low := mid + 1 end while_low_lt_high end for_c end for_b end for_a ;
found :
end
end.</lang>
- Output:
27^5 + 84^5 + 110^5 + 133^5 = 144^5
Arturo
<lang rebol>eulerSumOfPowers: function [top][
p5: map 0..top => [& ^ 5]
loop 4..top 'a [ loop 3..a-1 'b [ loop 2..b-1 'c [ loop 1..c-1 'd [ s: (get p5 a) + (get p5 b) + (get p5 c) + (get p5 d) if integer? index p5 s -> return ~"|a|^5 + |b|^5 + |c|^5 + |d|^5 = |index p5 s|^5" ] ] ] ]
return "not found" ; shouldn't reach here
]
print eulerSumOfPowers 249</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
AWK
<lang AWK>
- syntax: GAWK -f EULERS_SUM_OF_POWERS_CONJECTURE.AWK
BEGIN {
start_int = systime() main() printf("%d seconds\n",systime()-start_int) exit(0)
} function main( sum,s1,x0,x1,x2,x3) {
for (x0=1; x0<=250; x0++) { for (x1=1; x1<=x0; x1++) { for (x2=1; x2<=x1; x2++) { for (x3=1; x3<=x2; x3++) { sum = (x0^5) + (x1^5) + (x2^5) + (x3^5) s1 = int(sum ^ 0.2) if (sum == s1^5) { printf("%d^5 + %d^5 + %d^5 + %d^5 = %d^5\n",x0,x1,x2,x3,s1) return } } } } }
} </lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5 15 seconds
Bracmat
<lang bracmat> 0:?x0 & whl
' ( 1+!x0:<250:?x0 & out$(x0 !x0) & 0:?x1 & whl ' ( 1+!x1:~>!x0:?x1 & out$(x0 !x0 x1 !x1) & 0:?x2 & whl ' ( 1+!x2:~>!x1:?x2 & 0:?x3 & whl ' ( 1+!x3:~>!x2:?x3 & (!x0^5+!x1^5+!x2^5+!x3^5)^1/5 : ( #?y & out$(x0 !x0 x1 !x1 x2 !x2 x3 !x3 y !y) & 250:?x0:?x1:?x2:?x3 | ? ) ) ) ) )</lang>
Output
x0 133 x1 110 x2 84 x3 27 y 144
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
- include <stdio.h>
- 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#
Loops
<lang csharp>using System;
namespace EulerSumOfPowers {
class Program { const int MAX_NUMBER = 250;
static void Main(string[] args) { bool 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 (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 = Array.BinarySearch(fifth, sum); found = e >= 0; if (found) { Console.WriteLine("{0}^5 + {1}^5 + {2}^5 + {3}^5 = {4}^5", a + 1, b + 1, c + 1, d + 1, e + 1); } } } } } } }
}</lang>
Paired Powers, Mod 30, etc...
<lang csharp>using System; using System.Collections.Generic; using System.Linq; using System.Threading.Tasks;
namespace Euler_cs {
class Program { struct Pair { public int a, b; public Pair(int x, int y) { a = x; b = y; } }
static int min = 1, max = 250; static ulong[] p5; static SortedDictionary<ulong, Pair>[] sum2 = new SortedDictionary<ulong, Pair>[30];
static string Fmt(Pair p) { return string.Format("{0}^5 + {1}^5", p.a, p.b); }
public static void InitM() { for (int i = 0; i <= 29; i++) sum2[i] = new SortedDictionary<ulong, Pair>(); p5 = new ulong[max + 1]; p5[min] = Convert.ToUInt64(min) * Convert.ToUInt64(min); p5[min] *= p5[min] * Convert.ToUInt64(min); for (int i = min; i <= max - 1; i++) { for (int j = i + 1; j <= max; j++) { p5[j] = Convert.ToUInt64(j) * Convert.ToUInt64(j); p5[j] *= p5[j] * Convert.ToUInt64(j); if (j == max) continue; ulong x = p5[i] + p5[j]; sum2[x % 30].Add(x, new Pair(i, j)); } } }
static List<string> CalcM(int m) { List<string> res = new List<string>(); for (int i = max; i >= min; i--) { ulong p = p5[i]; int pm = i % 30, mp = (pm - m + 30) % 30; foreach (var s in sum2[m].Keys) { if (p <= s) break; ulong t = p - s; if (sum2[mp].Keys.Contains(t) && sum2[mp][t].a > sum2[m][s].b) res.Add(string.Format(" {1} + {2} = {0}^5", i, Fmt(sum2[m][s]), Fmt(sum2[mp][t]))); } } return res; }
static int Snip(string s) { int p = s.IndexOf("=") + 1; return Convert.ToInt32(s.Substring(p, s.IndexOf("^", p) - p)); }
static int CompareRes(string x, string y) { int res = Snip(x).CompareTo(Snip(y)); if (res == 0) res = x.CompareTo(y); return res; }
static int Validify(int def, string s) { int res = def, t = 0; int.TryParse(s, out t); if (t >= 1 && t < Math.Pow((double)(ulong.MaxValue >> 1), 0.2)) res = t; return res; }
static void Switch(ref int a, ref int b) { int t = a; a = b; b = t; }
static void Main(string[] args) { if (args.Count() > 1) { min = Validify(min, args[0]); max = Validify(max, args[1]); if (max < min) Switch(ref max, ref min); } else if (args.Count() == 1) max = Validify(max, args[0]); Console.WriteLine("Mod 30 shortcut with threading, checking from {0} to {1}...", min, max); List<string> res = new List<string>(); DateTime st = DateTime.Now; List<Task<List<string>>> taskList = new List<Task<List<string>>>(); InitM(); for (int j = 0; j <= 29; j++) { var jj = j; taskList.Add(Task.Run(() => CalcM(jj))); } Task.WhenAll(taskList); foreach (var item in taskList.Select(t => t.Result)) res.AddRange(item); res.Sort(CompareRes); foreach (var item in res) Console.WriteLine(item); Console.WriteLine(" Computation time to check entire space was {0} seconds", (DateTime.Now - st).TotalSeconds); if (System.Diagnostics.Debugger.IsAttached) Console.ReadKey(); } }
}</lang>
- Output:
(no command line arguments)
Mod 30 shortcut with threading, checking from 1 to 250...27^5 + 84^5 + 110^5 + 133^5 = 144^5
Computation time to check entire space was 0.0838058 seconds
(command line argument = "1000")
Mod 30 shortcut with threading, checking from 1 to 1000... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time to check entire space was 5.4109744 seconds
C++
First version
The simplest brute-force find is already reasonably quick: <lang cpp>#include <algorithm>
- include <iostream>
- include <cmath>
- include <set>
- include <vector>
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"; cout << "time=" << (int)((clock() - tm) * 1000 / CLOCKS_PER_SEC) << " milliseconds\r\n"; 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 cpp> 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. We also note that inside the innermost loop, we only need to search upward, so we hoist the downward linear search to the loop over x2. <lang cpp>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++) { 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) { cout << x0 << " " << x1 << " " << x2 << " " << x3 << " " << pow(sum, 1.0 / 5.0) << endl; return true; } } } } } // not found return false; }</lang> This reduces the timing to around 25 ms. We could equally well replace rs with an iterator inside pow5; the timing is unaffected.
For comparison with the C code, we also check the timing of an exhaustive search up to MAX=1000. (Don't try this in Python.) This takes 87.2 seconds on the same hardware, comparable to the results found by the C code authors, and supports their conclusion that the mod-30 trick used in the C solution leads to better scalability than the iterator optimizations.
Fortunately, we can incorporate the same trick into the above code, by inserting a forward jump to a feasible solution mod 30 into the loop over x3: <lang cpp> for (auto x3 = 1; x3 < x2; x3++) { // go straight to the first appropriate x3, mod 30 if (int err30 = (x0 + x1 + x2 + x3 - rs) % 30) x3 += 30 - err30; if (x3 >= x2) break; auto sum = s2 + pow5[x3];</lang> With this refinement, the exhaustive search up to MAX=1000 takes 16.9 seconds.
Thanks, C guys!
Second version
We can create a more efficient method by using the idea (taken from the EchoLisp listing below) of precomputing difference between pairs of fifth powers. If we combine this with the above idea of using linear search with memory, this still requires asymptotically O(N^4) time (because of the linear search within diffs), but is at least twice as fast as the solution above using the mod-30 trick. Exhaustive search up to MAX=1000 took 6.2 seconds for me (64-bit on 3.4GHz i7). It is not clear how it can be combined with the mod-30 trick.
The asymptotic behavior can be improved to O(N^3 ln N) by replacing the linear search with an increasing-increment "hunt" (and the outer linear search, which is also O(N^4), with a call to std::upper_bound). With this replacement, the first solution was found in 0.05 seconds; exhaustive search up to MAX=1000 took 2.80 seconds; and the second nontrivial solution (discarding multiples of the first solution), at y==2615, was found in 94.6 seconds. Note: there is no solution 2615, because 645^5 + 1523^5 + 1722^5 +2506^5 = 122 280 854 808 884 376, but 2615^5=122 280 854 808 884 375. This is an error due to limitation in mantissa of double type (52 bits). 128 bit type is required for the next solution 85359.
<lang cpp>template<class C_, class LT_> C_ Unique(const C_& src, const LT_& less)
{
C_ retval(src);
std::sort(retval.begin(), retval.end(), less);
retval.erase(unique(retval.begin(), retval.end()), retval.end());
return retval;
}
template<class I_, class P_> I_ HuntFwd(const I_& hint, const I_& end, const P_& less) // if less(x) is false, then less(x+1) must also be false { I_ retval(hint); int step = 1; // expanding phase while (end - retval > step) { I_ test = retval + step; if (!less(test)) break; retval = test; step <<= 1; } // contracting phase while (step > 1) { step >>= 1; if (end - retval <= step) continue; I_ test = retval + step; if (less(test)) retval = test; } if (retval != end && less(retval)) ++retval; return retval; }
bool DPFind(int how_many) { const int MAX = 1000; vector<double> pow5(MAX); for (int i = 1; i < MAX; i++) pow5[i] = (double)i * i * i * i * i; vector<pair<double, int>> diffs; for (int i = 2; i < MAX; ++i) { for (int j = 1; j < i; ++j) diffs.emplace_back(pow5[i] - pow5[j], j); } auto firstLess = [](const pair<double, int>& lhs, const pair<double, int>& rhs) { return lhs.first < rhs.first; }; diffs = Unique(diffs, firstLess);
for (int x4 = 4; x4 < MAX - 1; ++x4) { for (int x3 = 3; x3 < x4; ++x3) { // if (133 * x3 == 110 * x4) continue; // skip duplicates of first solution const auto s2 = pow5[x4] + pow5[x3]; auto pd = upper_bound(diffs.begin() + 1, diffs.end(), make_pair(s2, 0), firstLess) - 1; for (int x2 = 2; x2 < x3; ++x2) { const auto sum = s2 + pow5[x2]; pd = HuntFwd(pd, diffs.end(), [&](decltype(pd) it){ return it->first < sum; }); if (pd != diffs.end() && pd->first == sum && pd->second < x3) // find each solution only once { const double y = pow(pd->first + pow5[pd->second], 0.2); cout << x4 << " " << x3 << " " << x2 << " " << pd->second << " -> " << static_cast<int>(y + 0.5) << "\n"; if (--how_many <= 0) return true; } } } } return false; }</lang>
Thanks, EchoLisp guys!
Clojure
<lang Clojure> (ns test-p.core
(:require [clojure.math.numeric-tower :as math]) (:require [clojure.data.int-map :as i]))
(defn solve-power-sum [max-value max-sols]
" Finds solutions by using method approach of EchoLisp Large difference is we store a dictionary of all combinations of y^5 - x^5 with the x, y value so we can simply lookup rather than have to search " (let [pow5 (mapv #(math/expt % 5) (range 0 (* 4 max-value))) ; Pow5 = Generate Lookup table for x^5 y5-x3 (into (i/int-map) (for [x (range 1 max-value) ; For x0^5 + x1^5 + x2^5 + x3^5 = y^5 y (range (+ 1 x) (* 4 max-value))] ; compute y5-x3 = set of all possible differnences [(- (get pow5 y) (get pow5 x)) [x y]])) ; note: (get pow5 y) is closure for: pow5[y] solutions-found (atom 0)]
(for [x0 (range 1 max-value) ; Search over x0, x1, x2 for sums equal y5-x3 x1 (range 1 x0) x2 (range 1 x1) :when (< @solutions-found max-sols) :let [sum (apply + (map pow5 [x0 x1 x2]))] ; compute sum of items to the 5th power :when (contains? y5-x3 sum)] ; check if sum is in set of differences (do (swap! solutions-found inc) ; increment counter for solutions found (concat [x0 x1 x2] (get y5-x3 sum)))))) ; create result (since in set of differences)
- Output results with numbers in ascending order placing results into a set (i.e. #{}) so duplicates are discarded
; CPU i7 920 Quad Core @2.67 GHz clock Windows 10
(println (into #{} (map sort (solve-power-sum 250 1)))) ; MAX = 250, find only 1 value: Duration was 0.26 seconds (println (into #{} (map sort (solve-power-sum 1000 1000))));MAX = 1000, high max-value so all solutions found: Time = 4.8 seconds </lang> Output
1st Solution with MAX = 250 (Solution Time: 260 ms CPU i7 920 Quad Core) #{(27 84 110 133 144)) All Solutions with MAX = 1000 (Solution Time: 4.8 seconds CPU i7 920 Quad Core) #{(27 84 110 133 144) (162 504 660 798 864) (135 420 550 665 720) (108 336 440 532 576) (189 588 770 931 1008) (54 168 220 266 288) (81 252 330 399 432)}
COBOL
<lang cobol>
IDENTIFICATION DIVISION. PROGRAM-ID. EULER. DATA DIVISION. FILE SECTION. WORKING-STORAGE SECTION. 1 TABLE-LENGTH CONSTANT 250. 1 SEARCHING-FLAG PIC 9. 88 FINISHED-SEARCHING VALUE IS 1 WHEN SET TO FALSE IS 0. 1 CALC. 3 A PIC 999 USAGE COMPUTATIONAL-5. 3 B PIC 999 USAGE COMPUTATIONAL-5. 3 C PIC 999 USAGE COMPUTATIONAL-5. 3 D PIC 999 USAGE COMPUTATIONAL-5. 3 ABCD PIC 9(18) USAGE COMPUTATIONAL-5. 3 FIFTH-ROOT-OFFS PIC 999 USAGE COMPUTATIONAL-5. 3 POWER-COUNTER PIC 999 USAGE COMPUTATIONAL-5. 88 POWER-MAX VALUE TABLE-LENGTH. 1 PRETTY. 3 A PIC ZZ9. 3 FILLER VALUE "^5 + ". 3 B PIC ZZ9. 3 FILLER VALUE "^5 + ". 3 C PIC ZZ9. 3 FILLER VALUE "^5 + ". 3 D PIC ZZ9. 3 FILLER VALUE "^5 = ". 3 FIFTH-ROOT-OFFS PIC ZZ9. 3 FILLER VALUE "^5.".
1 FIFTH-POWER-TABLE OCCURS TABLE-LENGTH TIMES ASCENDING KEY IS FIFTH-POWER INDEXED BY POWER-INDEX. 3 FIFTH-POWER PIC 9(18) USAGE COMPUTATIONAL-5.
PROCEDURE DIVISION. MAIN-PARAGRAPH. SET FINISHED-SEARCHING TO FALSE. PERFORM POWERS-OF-FIVE-TABLE-INIT. PERFORM VARYING A IN CALC FROM 1 BY 1 UNTIL A IN CALC = TABLE-LENGTH
AFTER B IN CALC FROM 1 BY 1 UNTIL B IN CALC = A IN CALC
AFTER C IN CALC FROM 1 BY 1 UNTIL C IN CALC = B IN CALC
AFTER D IN CALC FROM 1 BY 1 UNTIL D IN CALC = C IN CALC
IF FINISHED-SEARCHING STOP RUN END-IF
PERFORM POWER-COMPUTATIONS
END-PERFORM.
POWER-COMPUTATIONS.
MOVE ZERO TO ABCD IN CALC.
ADD FIFTH-POWER(A IN CALC) FIFTH-POWER(B IN CALC) FIFTH-POWER(C IN CALC) FIFTH-POWER(D IN CALC) TO ABCD IN CALC.
SET POWER-INDEX TO 1.
SEARCH ALL FIFTH-POWER-TABLE WHEN FIFTH-POWER(POWER-INDEX) = ABCD IN CALC MOVE POWER-INDEX TO FIFTH-ROOT-OFFS IN CALC MOVE CORRESPONDING CALC TO PRETTY DISPLAY PRETTY END-DISPLAY SET FINISHED-SEARCHING TO TRUE END-SEARCH
EXIT PARAGRAPH.
POWERS-OF-FIVE-TABLE-INIT. PERFORM VARYING POWER-COUNTER FROM 1 BY 1 UNTIL POWER-MAX COMPUTE FIFTH-POWER(POWER-COUNTER) = POWER-COUNTER * POWER-COUNTER * POWER-COUNTER * POWER-COUNTER * POWER-COUNTER END-COMPUTE END-PERFORM. EXIT PARAGRAPH.
END PROGRAM EULER.
</lang> Output
133^5 + 110^5 + 84^5 + 27^5 = 144^5.
Common Lisp
<lang lisp> (ql:quickload :alexandria) (let ((fifth-powers (mapcar #'(lambda (x) (expt x 5))
(alexandria:iota 250)))) (loop named outer for x0 from 1 to (length fifth-powers) do (loop for x1 from 1 below x0 do (loop for x2 from 1 below x1 do (loop for x3 from 1 below x2 do (let ((x-sum (+ (nth x0 fifth-powers) (nth x1 fifth-powers) (nth x2 fifth-powers) (nth x3 fifth-powers)))) (if (member x-sum fifth-powers) (return-from outer (list x0 x1 x2 x3 (round (expt x-sum 0.2)))))))))))
</lang>
- Output:
(133 110 84 27 144)
D
First version
<lang d>import std.stdio, std.range, std.algorithm, std.typecons;
auto eulersSumOfPowers() {
enum maxN = 250; auto pow5 = iota(size_t(maxN)).map!(i => ulong(i) ^^ 5).array.assumeSorted;
foreach (immutable x0; 1 .. maxN) foreach (immutable x1; 1 .. x0) foreach (immutable x2; 1 .. x1) foreach (immutable x3; 1 .. x2) { immutable powSum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3]; if (pow5.contains(powSum)) return tuple(x0, x1, x2, x3, pow5.countUntil(powSum)); } assert(false);
}
void main() {
writefln("%d^5 + %d^5 + %d^5 + %d^5 == %d^5", eulersSumOfPowers[]);
}</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 == 144^5
Run-time about 0.64 seconds. A Range-based Haskell-like solution is missing because of Issue 14833.
Second version
<lang d>void main() {
import std.stdio, std.range, std.algorithm, std.typecons;
enum uint MAX = 250; uint[ulong] p5; Tuple!(uint, uint)[ulong] sum2;
foreach (immutable i; 1 .. MAX) { p5[ulong(i) ^^ 5] = i; foreach (immutable j; i .. MAX) sum2[ulong(i) ^^ 5 + ulong(j) ^^ 5] = tuple(i, j); }
const sk = sum2.keys.sort().release; foreach (p; p5.keys.sort()) foreach (immutable s; sk) { if (p <= s) break; if (p - s in sum2) { writeln(p5[p], " ", tuple(sum2[s][], sum2[p - s][])); return; // Finds first only. } }
}</lang>
- Output:
144 Tuple!(uint, uint, uint, uint)(27, 84, 110, 133)
Run-time about 0.10 seconds.
Third version
This solution is a brutal translation of the iterator-based C++ version, and it should be improved to use more idiomatic D Ranges. <lang d>import core.stdc.stdio, std.typecons, std.math, std.algorithm, std.range;
alias Pair = Tuple!(double, int); alias PairPtr = Pair*;
// If less(x) is false, then less(x + 1) must also be false. PairPtr huntForward(Pred)(PairPtr hint, const PairPtr end, const Pred less) pure nothrow @nogc {
PairPtr result = hint; int step = 1;
// Expanding phase. while (end - result > step) { PairPtr test = result + step; if (!less(test)) break; result = test; step <<= 1; }
// Contracting phase. while (step > 1) { step >>= 1; if (end - result <= step) continue; PairPtr test = result + step; if (less(test)) result = test; } if (result != end && less(result)) ++result; return result;
}
bool dPFind(int how_many) nothrow {
enum MAX = 1_000;
double[MAX] pow5; foreach (immutable i; 1 .. MAX) pow5[i] = double(i) ^^ 5;
Pair[] diffs0; // Will contain (MAX-1) * (MAX-2) / 2 pairs. foreach (immutable i; 2 .. MAX) foreach (immutable j; 1 .. i) diffs0 ~= Pair(pow5[i] - pow5[j], j);
// Remove pairs with duplicate first items. diffs0.length -= diffs0.sort!q{ a[0] < b[0] }.uniq.copy(diffs0).length; auto diffs = diffs0.assumeSorted!q{ a[0] < b[0] };
foreach (immutable x4; 4 .. MAX - 1) { foreach (immutable x3; 3 .. x4) { immutable s2 = pow5[x4] + pow5[x3]; auto pd0 = diffs[1 .. $].upperBound(Pair(s2, 0)); PairPtr pd = &pd0[0] - 1; foreach (immutable x2; 2 .. x3) { immutable sum = s2 + pow5[x2]; const PairPtr endPtr = &diffs[$ - 1] + 1; // This lambda heap-allocates. pd = huntForward(pd, endPtr, (in PairPtr p) pure => (*p)[0] < sum); if (pd != endPtr && (*pd)[0] == sum && (*pd)[1] < x3) { // Find each solution only once. immutable y = ((*pd)[0] + pow5[(*pd)[1]]) ^^ 0.2; printf("%d %d %d %d : %d\n", x4, x3, x2, (*pd)[1], cast(int)(y + 0.5)); if (--how_many <= 0) return true; } } } }
return false;
}
void main() nothrow {
if (!dPFind(100)) printf("Search finished.\n");
}</lang>
- Output:
133 110 27 84 : 144 133 110 84 27 : 144 266 220 54 168 : 288 266 220 168 54 : 288 399 330 81 252 : 432 399 330 252 81 : 432 532 440 108 336 : 576 532 440 336 108 : 576 665 550 135 420 : 720 665 550 420 135 : 720 798 660 162 504 : 864 798 660 504 162 : 864 Search finished.
Run-time about 7.1 seconds.
Delphi
See Pascal.
EasyLang
<lang>n = 250 len p5[] n len h5[] 65537 for i range n
p5[i] = i * i * i * i * i h5[p5[i] mod 65537] = 1
. func search a s . y .
y = -1 b = n while a + 1 < b i = (a + b) div 2 if p5[i] > s b = i elif p5[i] < s a = i else a = b y = i . .
. for x0 range n
for x1 range x0 sum1 = p5[x0] + p5[x1] for x2 range x1 sum2 = p5[x2] + sum1 for x3 range x2 sum = p5[x3] + sum2 if h5[sum mod 65537] = 1 call search x0 sum y if y >= 0 print x0 & " " & x1 & " " & x2 & " " & x3 & " " & y break 4 . . . . .
.</lang>
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>
Elixir
<lang elixir>defmodule Euler do
def sum_of_power(max \\ 250) do {p5, sum2} = setup(max) sk = Enum.sort(Map.keys(sum2)) Enum.reduce(Enum.sort(Map.keys(p5)), Map.new, fn p,map -> sum(sk, p5, sum2, p, map) end) end defp setup(max) do Enum.reduce(1..max, {%{}, %{}}, fn i,{p5,sum2} -> i5 = i*i*i*i*i add = for j <- i..max, into: sum2, do: {i5 + j*j*j*j*j, [i,j]} {Map.put(p5, i5, i), add} end) end defp sum([], _, _, _, map), do: map defp sum([s|_], _, _, p, map) when p<=s, do: map defp sum([s|t], p5, sum2, p, map) do if sum2[p - s], do: sum(t, p5, sum2, p, Map.put(map, Enum.sort(sum2[s] ++ sum2[p-s]), p5[p])), else: sum(t, p5, sum2, p, map) end
end
Enum.each(Euler.sum_of_power, fn {k,v} ->
IO.puts Enum.map_join(k, " + ", fn i -> "#{i}**5" end) <> " = #{v}**5"
end)</lang>
- Output:
27**5 + 84**5 + 110**5 + 133**5 = 144**5
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"
Factor
This solution uses Factor's backtrack
vocabulary (based on continuations) to simplify the reduction of the search space. Each time xn
is called, a new summand is introduced which can only take on a value as high as the previous summand - 1. This also creates a checkpoint for the backtracker. fail
causes the backtracking to occur.
<lang factor>USING: arrays backtrack kernel literals math.functions
math.ranges prettyprint sequences ;
CONSTANT: pow5 $[ 0 250 [a,b) [ 5 ^ ] map ]
- xn ( n1 -- n2 n2 ) [1,b) amb-lazy dup ;
250 xn xn xn xn drop 4array dup pow5 nths sum dup pow5 member? [ pow5 index suffix . ] [ 2drop fail ] if</lang>
- Output:
{ 133 110 84 27 144 }
Forth
<lang forth>
- sq dup * ;
- 5^ dup sq sq * ;
create pow5 250 cells allot
- noname
250 0 DO i 5^ pow5 i cells + ! LOOP ; execute
- @5^ cells pow5 + @ ;
- solution? ( n -- n )
pow5 250 cells bounds DO dup i @ = IF drop i pow5 - cell / unloop EXIT THEN cell +LOOP drop 0 ;
\ GFORTH only provides 2 index variables: i, j \ so the code creates locals for two outer loop vars, k & l
- euler ( -- )
250 4 DO i { l } l 3 DO i { k } k 2 DO i 1 DO i @5^ j @5^ + k @5^ + l @5^ + solution? dup IF l . k . j . i . . cr unloop unloop unloop unloop EXIT ELSE drop THEN LOOP LOOP LOOP LOOP ;
euler bye </lang>
- Output:
$ gforth-fast ./euler.fs 133 110 84 27 144
Fortran
FORTRAN IV
To solve this problem, we must handle integers up 250**5 ~= 9.8*10**11 . So we need integers with at less 41 bits. In 1966 all Fortrans were not equal. On IBM360, INTEGER was a 32-bit integer; on CDC6600, INTEGER was a 60-bit integer. And Leon J. Lander and Thomas R. Parkin used the CDC6600. <lang fortran>C EULER SUM OF POWERS CONJECTURE - FORTRAN IV C FIND I1,I2,I3,I4,I5 : I1**5+I2**5+I3**5+I4**5=I5**5
INTEGER I,P5(250),SUMX MAXN=250 DO 1 I=1,MAXN 1 P5(I)=I**5 DO 6 I1=1,MAXN DO 6 I2=1,MAXN DO 6 I3=1,MAXN DO 6 I4=1,MAXN SUMX=P5(I1)+P5(I2)+P5(I3)+P5(I4) I5=1 2 IF(I5-MAXN) 3,3,6 3 IF(P5(I5)-SUMX) 5,4,6 4 WRITE(*,300) I1,I2,I3,I4,I5 STOP 5 I5=I5+1 GOTO 2 6 CONTINUE 300 FORMAT(5(1X,I3)) END </lang>
- Output:
27 84 110 133 144
Fortran 95
<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
FreeBASIC
<lang freebasic>' version 14-09-2015 ' compile with: fbc -s console
' some constants calculated when the program is compiled
Const As UInteger max = 250 Const As ULongInt pow5_max = CULngInt(max) * max * max * max * max ' limit x1, x2, x3 Const As UInteger limit_x1 = (pow5_max / 4) ^ 0.2 Const As UInteger limit_x2 = (pow5_max / 3) ^ 0.2 Const As UInteger limit_x3 = (pow5_max / 2) ^ 0.2
' ------=< MAIN >=------
Dim As ULongInt pow5(max), ans1, ans2, ans3 Dim As UInteger x1, x2, x3, x4, x5 , m1, m2
Cls : Print
For x1 = 1 To max
pow5(x1) = CULngInt(x1) * x1 * x1 * x1 * x1
Next x1
For x1 = 1 To limit_x1
For x2 = x1 +1 To limit_x2 m1 = x1 + x2 ans1 = pow5(x1) + pow5(x2) If ans1 > pow5_max Then Exit For For x3 = x2 +1 To limit_x3 ans2 = ans1 + pow5(x3) If ans2 > pow5_max Then Exit For m2 = (m1 + x3) Mod 30 If m2 = 0 Then m2 = 30 For x4 = x3 +1 To max -1 ans3 = ans2 + pow5(x4) If ans3 > pow5_max Then Exit For For x5 = x4 + m2 To max Step 30 If ans3 < pow5(x5) Then Exit For If ans3 = pow5(x5) Then Print x1; "^5 + "; x2; "^5 + "; x3; "^5 + "; _ x4; "^5 = "; x5; "^5" Exit For, For EndIf Next x5 Next x4 Next x3 Next x2
Next x1
Print Print "done"
' empty keyboard buffer While Inkey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
27^5 + 84^5 + 110^5 + 133^5 = 144^5
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
Groovy
<lang groovy>class EulerSumOfPowers {
static final int MAX_NUMBER = 250
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 (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 = Arrays.binarySearch(fifth, sum) found = (e >= 0) if (found) { println("${a + 1}^5 + ${b + 1}^5 + ${c + 1}^5 + ${d + 1}^5 + ${e + 1}^5") } } } } } }
}</lang>
- Output:
27^5 + 84^5 + 110^5 + 133^5 + 144^5
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)
Or, using dictionaries of powers and sums, and thus rather faster:
<lang haskell>import qualified Data.Map.Strict as M import Data.List (find, intercalate) import Data.Maybe (maybe)
EULER'S SUM OF POWERS CONJECTURE -----------
counterExample :: (M.Map Int (Int, Int), M.Map Int Int) -> Maybe (Int, Int) counterExample (sumMap, powerMap) =
find (\(p, s) -> M.member (p - s) sumMap) (M.keys powerMap >>= (((>>=) . flip takeWhile (M.keys sumMap) . (>)) <*> \ p s -> [(p, s)]))
sumMapForRange :: [Int] -> M.Map Int (Int, Int) sumMapForRange xs =
M.fromList [ ((x ^ 5) + (y ^ 5), (x, y)) | x <- xs , y <- tail xs , x > y ]
powerMapForRange :: [Int] -> M.Map Int Int powerMapForRange = M.fromList . (zip =<< fmap (^ 5))
TEST -------------------------
main :: IO () main =
putStrLn $ "Euler's sum of powers conjecture – " <> maybe ("no counter-example found in the range " <> rangeString xs) (showExample sumsAndPowers xs) (counterExample sumsAndPowers) where xs = [1 .. 249] sumsAndPowers = ((,) . sumMapForRange <*> powerMapForRange) xs
showExample :: (M.Map Int (Int, Int), M.Map Int Int) -> [Int] -> (Int, Int) -> String showExample (sumMap, powerMap) xs (p, s) =
"a counter-example in range " <> rangeString xs <> ":\n\n" <> intercalate "^5 + " (show <$> [a, b, c, d]) <> "^5 = " <> show (powerMap M.! p) <> "^5" where (a, b) = sumMap M.! (p - s) (c, d) = sumMap M.! s
rangeString :: [Int] -> String rangeString [] = "[]" rangeString (x:xs) = '[' : show x <> " .. " <> show (last xs) <> "]"</lang>
- Output:
Euler's sum of powers conjecture – a counter-example in range [1 .. 249]: 133^5 + 110^5 + 84^5 + 27^5 = 144^5
J
<lang J> require 'stats'
(#~ (= <.)@((+/"1)&.:(^&5)))1+4 comb 248
27 84 110 133</lang>
Explanation:
<lang J>1+4 comb 248</lang> finds all the possibilities for our four arguments.
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.
Here's a significantly faster approach (about 100 times faster), based on the echolisp implementation:
<lang J>find5=:3 :0
y=. 250 n=. i.y p=. n^5 a=. (#~ 0&<),-/~p s=. /:~a l=. (i.*:y)(#~ 0&<),-/~p c=. 3 comb <.5%:(y^5)%4 t=. +/"1 c{p x=. (t e. s)#t |.,&<&~./|:(y,y)#:l#~a e. x
)</lang>
Use:
<lang J> find5 ┌─────────────┬───┐ │27 84 110 133│144│ └─────────────┴───┘</lang>
Note that this particular implementation is a bit hackish, since it relies on the solution being unique for the range of numbers being considered. If there were more solutions it would take a little extra code (though not much time) to untangle them.
Java
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
ES5
<lang javascript>var eulers_sum_of_powers = function (iMaxN) {
var aPow5 = []; var oPow5ToN = {};
for (var iP = 0; 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,l 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>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
This
that verify: a^5 + b^5 + c^5 + d^5 = x^5
<lang javascript>var N=1000, first=false
var ns={}, npv=[]
for (var n=0; n<=N; n++) {
var np=Math.pow(n,5); ns[np]=n; npv.push(np)
}
loop:
for (var a=1; a<=N; a+=1)
for (var b=a+1; b<=N; b+=1)
for (var c=b+1; c<=N; c+=1)
for (var d=c+1; d<=N; d+=1) {
var x = ns[ npv[a]+npv[b]+npv[c]+npv[d] ]
if (!x) continue
print( [a, b, c, d, x] )
if (first) break loop
}
function print(c) {
var e='5', ep=e+' + '
document.write(c[0], ep, c[1], ep, c[2], ep, c[3], e, ' = ', c[4], e, '
')
}</lang>
Or this
that verify: a^5 + b^5 + c^5 + d^5 = x^5
<lang javascript>var N=1000, first=false var npv=[], M=30 // x^5 == x modulo M (=2*3*5) for (var n=0; n<=N; n+=1) npv[n]=Math.pow(n, 5) var mx=1+npv[N]; while(n<=N+M) npv[n++]=mx
loop:
for (var a=1; a<=N; a+=1)
for (var b=a+1; b<=N; b+=1)
for (var c=b+1; c<=N; c+=1)
for (var t=npv[a]+npv[b]+npv[c], d=c+1, x=t%M+d; (n=t+npv[d])<mx; d+=1, x+=1) {
while (npv[x]<=n) x+=M; x-=M // jump over M=30 values for x>d
if (npv[x] != n) continue
print( [a, b, c, d, x] )
if (first) break loop;
}
function print(c) {
var e='5', ep=e+' + '
document.write(c[0], ep, c[1], ep, c[2], ep, c[3], e, ' = ', c[4], e, '
')
}</lang>
Or this
that verify: a^5 + b^5 + c^5 = x^5 - d^5
<lang javascript>var N=1000, first=false
var dxs={}, pow=Math.pow
for (var d=1; d<=N; d+=1)
for (var dp=pow(d,5), x=d+1; x<=N; x+=1)
dxs[pow(x,5)-dp]=[d,x]
loop:
for (var a=1; a<N; a+=1)
for (var ap=pow(a,5), b=a+1; b<N; b+=1)
for (var abp=ap+pow(b,5), c=b+1; c<N; c+=1) {
var dx = dxs[ abp+pow(c,5) ]
if (!dx || c >= dx[0]) continue
print( [a, b, c].concat( dx ) )
if (first) break loop
}
function print(c) {
var e='5', ep=e+' + '
document.write(c[0], ep, c[1], ep, c[2], ep, c[3], e, ' = ', c[4], e, '
')
}</lang>
Or this
that verify: a^5 + b^5 = x^5 - (c^5 + d^5)
<lang javascript>var N=1000, first=false
var is={}, ipv=[], ijs={}, ijpv=[], pow=Math.pow
for (var i=1; i<=N; i+=1) {
var ip=pow(i,5); is[ip]=i; ipv.push(ip)
for (var j=i+1; j<=N; j+=1) {
var ijp=ip+pow(j,5); ijs[ijp]=[i,j]; ijpv.push(ijp)
}
}
ijpv.sort( function (a,b) {return a - b } )
loop:
for (var i=0, ei=ipv.length; i<ei; i+=1)
for (var xp=ipv[i], j=0, je=ijpv.length; j<je; j+=1) {
var cdp = ijpv[j]
if (cdp >= xp) break
var cd = ijs[xp-cdp]
if (!cd) continue
var ab = ijs[cdp]
if (ab[1] >= cd[0]) continue
print( [].concat(ab, cd, is[xp]) )
if (first) break loop
}
function print(c) {
var e='5', ep=e+' + '
document.write(c[0], ep, c[1], ep, c[2], ep, c[3], e, ' = ', c[4], e, '
')
}</lang>
- Output:
275 + 845 + 1105 + 1335 = 1445 545 + 1685 + 2205 + 2665 = 2885 815 + 2525 + 3305 + 3995 = 4325 1085 + 3365 + 4405 + 5325 = 5765 1355 + 4205 + 5505 + 6655 = 7205 1625 + 5045 + 6605 + 7985 = 8645
ES6
Procedural
<lang JavaScript>(() => {
'use strict';
const eulersSumOfPowers = intMax => { const pow = Math.pow, xs = range(0, intMax) .map(x => pow(x, 5)), dct = xs.reduce((a, x, i) => (a[x] = i, a ), {});
for (let a = 1; a <= intMax; a++) { for (let b = 2; b <= a; b++) { for (let c = 3; c <= b; c++) { for (let d = 4; d <= c; d++) { const sumOfPower = dct[xs[a] + xs[b] + xs[c] + xs[d]]; if (sumOfPower !== undefined) { return [a, b, c, d, sumOfPower]; } } } } } return undefined; };
// range :: Int -> Int -> [Int] const range = (m, n) => Array.from({ length: Math.floor(n - m) + 1 }, (_, i) => m + i);
// TEST const soln = eulersSumOfPowers(250); return soln ? soln.slice(0, 4) .map(x => `${x}^5`) .join(' + ') + ` = ${soln[4]}^5` : 'No solution found.'
})();</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
Functional
Using dictionaries of powers and sums, and a little faster than the procedural version above:
<lang javascript>(() => {
'use strict';
const main = () => {
const iFrom = 1, iTo = 249, xs = enumFromTo(1, 249), p5 = x => Math.pow(x, 5);
const // powerMap :: Dict Int Int powerMap = mapFromList( zip(map(p5, xs), xs) ), // sumMap :: Dict Int (Int, Int) sumMap = mapFromList( bind( xs, x => bind( tail(xs), y => Tuple( p5(x) + p5(y), Tuple(x, y) ) ) ) );
// mbExample :: Maybe (Int, Int) const mbExample = find( tpl => member(fst(tpl) - snd(tpl), sumMap), bind( map(x => parseInt(x, 10), keys(powerMap) ), p => bind( takeWhile( x => x < p, map(x => parseInt(x, 10), keys(sumMap) ) ), s => [Tuple(p, s)] ) ) );
// showExample :: (Int, Int) -> String const showExample = tpl => { const [p, s] = Array.from(tpl); const [a, b] = Array.from(sumMap[p - s]); const [c, d] = Array.from(sumMap[s]); return 'Counter-example found:\n' + intercalate( '^5 + ', map(str, [a, b, c, d]) ) + '^5 = ' + str(powerMap[p]) + '^5'; };
return maybe( 'No counter-example found', showExample, mbExample ); };
// GENERIC FUNCTIONS ----------------------------------
// Just :: a -> Maybe a const Just = x => ({ type: 'Maybe', Nothing: false, Just: x });
// Nothing :: Maybe a const Nothing = () => ({ type: 'Maybe', Nothing: true, });
// Tuple (,) :: a -> b -> (a, b) const Tuple = (a, b) => ({ type: 'Tuple', '0': a, '1': b, length: 2 });
// bind (>>=) :: [a] -> (a -> [b]) -> [b] const bind = (xs, mf) => [].concat.apply([], xs.map(mf));
// concat :: a -> [a] // concat :: [String] -> String const concat = xs => 0 < xs.length ? (() => { const unit = 'string' !== typeof xs[0] ? ( [] ) : ; return unit.concat.apply(unit, xs); })() : [];
// enumFromTo :: (Int, Int) -> [Int] const enumFromTo = (m, n) => Array.from({ length: 1 + n - m }, (_, i) => m + i);
// find :: (a -> Bool) -> [a] -> Maybe a const find = (p, xs) => { for (let i = 0, lng = xs.length; i < lng; i++) { if (p(xs[i])) return Just(xs[i]); } return Nothing(); };
// fst :: (a, b) -> a const fst = tpl => tpl[0];
// intercalate :: [a] -> a -> [a] // intercalate :: String -> [String] -> String const intercalate = (sep, xs) => 0 < xs.length && 'string' === typeof sep && 'string' === typeof xs[0] ? ( xs.join(sep) ) : concat(intersperse(sep, xs));
// intersperse(0, [1,2,3]) -> [1, 0, 2, 0, 3]
// intersperse :: a -> [a] -> [a] // intersperse :: Char -> String -> String const intersperse = (sep, xs) => { const bln = 'string' === typeof xs; return xs.length > 1 ? ( (bln ? concat : x => x)( (bln ? ( xs.split() ) : xs) .slice(1) .reduce((a, x) => a.concat([sep, x]), [xs[0]]) )) : xs; };
// keys :: Dict -> [String] const keys = Object.keys;
// Returns Infinity over objects without finite length. // This enables zip and zipWith to choose the shorter // argument when one is non-finite, like cycle, repeat etc
// length :: [a] -> Int const length = xs => (Array.isArray(xs) || 'string' === typeof xs) ? ( xs.length ) : Infinity;
// map :: (a -> b) -> [a] -> [b] const map = (f, xs) => (Array.isArray(xs) ? ( xs ) : xs.split()).map(f);
// mapFromList :: [(k, v)] -> Dict const mapFromList = kvs => kvs.reduce( (a, kv) => { const k = kv[0]; return Object.assign(a, { [ (('string' === typeof k) && k) || JSON.stringify(k) ]: kv[1] }); }, {} );
// Default value (v) if m.Nothing, or f(m.Just)
// maybe :: b -> (a -> b) -> Maybe a -> b const maybe = (v, f, m) => m.Nothing ? v : f(m.Just);
// member :: Key -> Dict -> Bool const member = (k, dct) => k in dct;
// snd :: (a, b) -> b const snd = tpl => tpl[1];
// str :: a -> String const str = x => x.toString();
// tail :: [a] -> [a] const tail = xs => 0 < xs.length ? xs.slice(1) : [];
// take :: Int -> [a] -> [a] // take :: Int -> String -> String const take = (n, xs) => 'GeneratorFunction' !== xs.constructor.constructor.name ? ( xs.slice(0, n) ) : [].concat.apply([], Array.from({ length: n }, () => { const x = xs.next(); return x.done ? [] : [x.value]; }));
// takeWhile :: (a -> Bool) -> [a] -> [a] // takeWhile :: (Char -> Bool) -> String -> String const takeWhile = (p, xs) => xs.constructor.constructor.name !== 'GeneratorFunction' ? (() => { const lng = xs.length; return 0 < lng ? xs.slice( 0, until( i => lng === i || !p(xs[i]), i => 1 + i, 0 ) ) : []; })() : takeWhileGen(p, xs);
// until :: (a -> Bool) -> (a -> a) -> a -> a const until = (p, f, x) => { let v = x; while (!p(v)) v = f(v); return v; };
// Use of `take` and `length` here allows for zipping with non-finite // lists - i.e. generators like cycle, repeat, iterate.
// zip :: [a] -> [b] -> [(a, b)] const zip = (xs, ys) => { const lng = Math.min(length(xs), length(ys)); return Infinity !== lng ? (() => { const bs = take(lng, ys); return take(lng, xs).map((x, i) => Tuple(x, bs[i])); })() : zipGen(xs, ys); };
// MAIN --- return main();
})();</lang>
- Output:
Counter-example found: 133^5 + 110^5 + 84^5 + 27^5 = 144^5
jq
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
- 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>
The task: <lang jq>sum_of_powers_conjecture(249)</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.
Kotlin
<lang scala>fun main(args: Array<String>) {
val p5 = LongArray(250){ it.toLong() * it * it * it * it } var sum: Long var y: Int var found = false loop@ for (x0 in 0 .. 249) for (x1 in 0 .. x0 - 1) for (x2 in 0 .. x1 - 1) for (x3 in 0 .. x2 - 1) { sum = p5[x0] + p5[x1] + p5[x2] + p5[x3] y = p5.binarySearch(sum) if (y >= 0) { println("$x0^5 + $x1^5 + $x2^5 + $x3^5 = $y^5") found = true break@loop } } if (!found) println("No solution was found")
}</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
Lua
Brute force but still takes under two seconds with LuaJIT. <lang Lua>-- Fast table search (only works if table values are in order) function binarySearch (t, n)
local start, stop, mid = 1, #t while start < stop do mid = math.floor((start + stop) / 2) if n == t[mid] then return mid elseif n < t[mid] then stop = mid - 1 else start = mid + 1 end end return nil
end
-- Test Euler's sum of powers conjecture function euler (limit)
local pow5, sum = {} for i = 1, limit do pow5[i] = i^5 end for x0 = 1, limit do for x1 = 1, x0 do for x2 = 1, x1 do for x3 = 1, x2 do sum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3] if binarySearch(pow5, sum) then print(x0 .. "^5 + " .. x1 .. "^5 + " .. x2 .. "^5 + " .. x3 .. "^5 = " .. sum^(1/5) .. "^5") return true end end end end end return false
end
-- Main procedure if euler(249) then
print("Time taken: " .. os.clock() .. " seconds")
else
print("Looks like he was right after all...")
end</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5 Time taken: 1.247 seconds
Mathematica /Wolfram Language
<lang Mathematica>Sort[FindInstance[
x0^5 + x1^5 + x2^5 + x3^5 == y^5 && x0 > 0 && x1 > 0 && x2 > 0 && x3 > 0, {x0, x1, x2, x3, y}, Integers]1, All, -1]</lang>
- Output:
{27,84,110,133,144}
Microsoft Small Basic
<lang smallbasic>' Euler sum of powers conjecture - 03/07/2015
'find: x1^5+x2^5+x3^5+x4^5=x5^5 '-> x1=27 x2=84 x3=110 x4=133 x5=144 maxn=250 For i=1 to maxn
p5[i]=Math.Power(i,5)
EndFor For x1=1 to maxn-4 For x2=x1+1 to maxn-3 'TextWindow.WriteLine("x1="+x1+", x2="+x2) For x3=x2+1 to maxn-2 'TextWindow.WriteLine("x1="+x1+", x2="+x2+", x3="+x3) For x4=x3+1 to maxn-1 'TextWindow.WriteLine("x1="+x1+", x2="+x2+", x3="+x3+", x4="+x4) x5=x4+1 valx=p5[x5] sumx=p5[x1]+p5[x2]+p5[x3]+p5[x4] While x5<=maxn and valx<=sumx If valx=sumx Then TextWindow.WriteLine("Found!") TextWindow.WriteLine("-> "+x1+"^5+"+x2+"^5+"+x3+"^5+"+x4+"^5="+x5+"^5") TextWindow.WriteLine("x5^5="+sumx) Goto EndPgrm EndIf x5=x5+1 valx=p5[x5] EndWhile 'x5 EndFor 'x4 EndFor 'x3 EndFor 'x2 EndFor 'x1 EndPgrm: </lang>
- Output:
Found! -> 27^5+84^5+110^5+133^5=144^5 x5^5=61917364224
Modula-2
<lang modula2>MODULE EulerConjecture;
FROM FormatString IMPORT FormatString;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
PROCEDURE Pow5(a : LONGINT) : LONGINT; BEGIN
RETURN a * a * a * a * a
END Pow5;
VAR
buf : ARRAY[0..63] OF CHAR; a,b,c,d,e,sum,curr : LONGINT;
BEGIN
FOR a:=0 TO 250 DO FOR b:=a TO 250 DO IF b=a THEN CONTINUE END; FOR c:=b TO 250 DO IF (c=a) OR (c=b) THEN CONTINUE END; FOR d:=c TO 250 DO IF (d=a) OR (d=b) OR (d=c) THEN CONTINUE END; sum := Pow5(a) + Pow5(b) + Pow5(c) + Pow5(d); FOR e:=d TO 250 DO IF (e=a) OR (e=b) OR (e=c) OR (e=d) THEN CONTINUE END; curr := Pow5(e); IF (sum#0) AND (sum=curr) THEN FormatString("%l^5 + %l^5 + %l^5 + %l^5 = %l^5\n", buf, a, b, c, d, e); WriteString(buf) ELSIF curr > sum THEN BREAK END END; END; END; END; END;
WriteString("Done"); WriteLn; ReadChar
END EulerConjecture.</lang>
Nim
<lang Nim>
- Brute force approach
import times
- assumes an array of non-decreasing positive integers
proc binarySearch(a : openArray[int], target : int) : int =
var left, right, mid : int left = 0 right = len(a) - 1 while true : if left > right : return 0 # no match found mid = (left + right) div 2 if a[mid] < target : left = mid + 1 elif a[mid] > target : right = mid - 1 else : return mid # match found
var
p5 : array[250, int] sum = 0 y, t1 : int
let t0 = cpuTime()
for i in 1 .. 249 :
p5[i] = i * i * i * i * i
for x0 in 1 .. 249 :
for x1 in 1 .. x0 - 1 : for x2 in 1 .. x1 - 1 : for x3 in 1 .. x2 - 1 : sum = p5[x0] + p5[x1] + p5[x2] + p5[x3] y = binarySearch(p5, sum) if y > 0 : t1 = int((cputime() - t0) * 1000.0) echo "Time : ", t1, " milliseconds" echo $x0 & "^5 + " & $x1 & "^5 + " & $x2 & "^5 + " & $x3 & "^5 = " & $y & "^5" quit()
if y == 0 :
echo "No solution was found"
</lang>
- Output:
Time : 156 milliseconds 133^5 + 110^5 + 84^5 + 27^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]
PARI/GP
Naive script: <lang parigp>forvec(v=vector(4,i,[0,250]), if(ispower(v[1]^5+v[2]^5+v[3]^5+v[4]^5,5,&n), print(n" "v)), 2)</lang>
- Output:
144 [27, 84, 110, 133]
Naive + caching (setbinop
):
<lang parigp>{
v2=setbinop((x,y)->[min(x,y),max(x,y),x^5+y^5],[0..250]); \\ sums of two fifth powers
for(i=2,#v2,
for(j=1,i-1, if(v2[i][2]<v2[j][2] && ispower(v2[i][3]+v2[j][3],5,&n) && #(v=Set([v2[i][1],v2[i][2],v2[j][1],v2[j][2]]))==4, print(n" "v) ) )
) }</lang>
- Output:
144 [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
Brute Force: <lang perl>use constant MAX => 250; my @p5 = (0,map { $_**5 } 1 .. MAX-1); my $s = 0; my %p5 = map { $_ => $s++ } @p5; for my $x0 (1..MAX-1) {
for my $x1 (1..$x0-1) { for my $x2 (1..$x1-1) { for my $x3 (1..$x2-1) { my $sum = $p5[$x0] + $p5[$x1] + $p5[$x2] + $p5[$x3]; die "$x3 $x2 $x1 $x0 $p5{$sum}\n" if exists $p5{$sum}; } } }
}</lang>
- Output:
27 84 110 133 144
Adding some optimizations makes it 5x faster with similar output, but obfuscates things.
<lang perl>use constant MAX => 250;
my @p5 = (0,map { $_**5 } 1 .. MAX-1); my $rs = 5; for my $x0 (1..MAX-1) {
for my $x1 (1..$x0-1) { for my $x2 (1..$x1-1) { my $s2 = $p5[$x0] + $p5[$x1] + $p5[$x2]; $rs-- while $rs > 0 && $p5[$rs] > $s2; for (my $x3 = 1; $x3 < $x2; $x3++) { my $e30 = ($x0 + $x1 + $x2 + $x3 - $rs) % 30; $x3 += (30-$e30) if $e30; last if $x3 >= $x2; my $sum = $s2 + $p5[$x3]; $rs++ while $rs < MAX-1 && $p5[$rs] < $sum; die "$x3 $x2 $x1 $x0 $rs\n" if $p5[$rs] == $sum; } } }
}</lang>
Phix
Around four seconds, not spectacularly fast. My naive brute force was over a minute. This is not where Phix shines.
Quitting when the first is found drops the main loop to 0.7s, so 1.1s in all, vs 4.3s for the full search.
Without the return 0, you just get six permutes (of ordered pairs) for 144.
<lang Phix>constant MAX = 250
constant p5 = new_dict(),
sum2 = new_dict()
atom t0 = time() for i=1 to MAX do
atom i5 = power(i,5) setd(i5,i,p5) for j=1 to i-1 do atom j5 = power(j,5) setd(j5+i5,{j,i},sum2) end for
end for
?time()-t0
function forsum2(object s, object data, object p)
if p<=s then return 0 end if integer k = getd_index(p-s,sum2) if k!=NULL then ?getd(p,p5)&data&getd_by_index(k,sum2) return 0 -- (show one solution per p) end if return 1
end function
function forp5(object key, object /*data*/, object /*user_data*/)
traverse_dict(routine_id("forsum2"),key,sum2) return 1
end function
traverse_dict(routine_id("forp5"),0,p5)
?time()-t0</lang>
- Output:
0.421 {144,27,84,110,133} 4.312
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
PicoLisp
<lang PicoLisp>(off P) (off S)
(for I 250
(idx 'P (list (setq @@ (** I 5)) I) T ) (for (J I (>= 250 J) (inc J)) (idx 'S (list (+ @@ (** J 5)) (list I J)) T ) ) )
(println
(catch 'found (for A (idx 'P) (for B (idx 'S) (T (<= (car A) (car B))) (and (lup S (- (car A) (car B))) (throw 'found (conc (cadr (lup S (car B))) (cadr (lup S (- (car A) (car B)))) (cdr (lup P (car A))) ) ) ) ) ) ) )</lang>
- Output:
(27 84 110 133 144)
PowerShell
Brute Force Search
This is a slow algorithm, so attempts have been made to speed it up, including pre-computing the powers, using an ArrayList for them, and using [int] to cast the 5th root rather than use truncate.
<lang powershell># EULER.PS1
$max = 250
$powers = New-Object System.Collections.ArrayList for ($i = 0; $i -lt $max; $i++) {
$tmp = $powers.Add([Math]::Pow($i, 5))
}
for ($x0 = 1; $x0 -lt $max; $x0++) {
for ($x1 = 1; $x1 -lt $x0; $x1++) { for ($x2 = 1; $x2 -lt $x1; $x2++) { for ($x3 = 1; $x3 -lt $x2; $x3++) { $sum = $powers[$x0] + $powers[$x1] + $powers[$x2] + $powers[$x3] $S1 = [int][Math]::pow($sum,0.2)
if ($sum -eq $powers[$S1]) { Write-host "$x0^5 + $x1^5 + $x2^5 + $x3^5 = $S1^5" return } } } }
}</lang>
- Output:
PS > measure-command { .\euler.ps1 | out-default } 133^5 + 110^5 + 84^5 + 27^5 = 144^5 Days : 0 Hours : 0 Minutes : 0 Seconds : 31 Milliseconds : 608 Ticks : 316082251 TotalDays : 0.000365835938657407
Prolog
<lang prolog> makepowers :-
retractall(pow5(_, _)), between(1, 249, X), Y is X * X * X * X * X, assert(pow5(X, Y)), fail.
makepowers.
within(A, Bx, N) :- % like between but with an exclusive upper bound
succ(B, Bx), between(A, B, N).
solution(X0, X1, X2, X3, Y) :-
makepowers, within(4, 250, X0), pow5(X0, X0_5th), within(3, X0, X1), pow5(X1, X1_5th), within(2, X1, X2), pow5(X2, X2_5th), within(1, X2, X3), pow5(X3, X3_5th), Y_5th is X0_5th + X1_5th + X2_5th + X3_5th, pow5(Y, Y_5th).
</lang>
- Output:
?- solution(X0,X1,X2,X3,Y). X0 = 133, X1 = 110, X2 = 84, X3 = 27, Y = 144 .
PureBasic
<lang PureBasic> EnableExplicit
- assumes an array of non-decreasing positive integers
Procedure.q BinarySearch(Array a.q(1), Target.q)
Protected l = 0, r = ArraySize(a()), m Repeat If l > r : ProcedureReturn 0 : EndIf; no match found m = (l + r) / 2 If a(m) < target l = m + 1 ElseIf a(m) > target r = m - 1 Else ProcedureReturn m ; match found EndIf ForEver
EndProcedure
Define i, x0, x1, x2, x3, y Define.q sum Define Dim p5.q(249)
For i = 1 To 249
p5(i) = i * i * i * i * i
Next
If OpenConsole()
For x0 = 1 To 249 For x1 = 1 To x0 - 1 For x2 = 1 To x1 - 1 For x3 = 1 To x2 - 1 sum = p5(x0) + p5(x1) + p5(x2) + p5(x3) y = BinarySearch(p5(), sum) If y > 0 PrintN(Str(x0) + "^5 + " + Str(x1) + "^5 + " + Str(x2) + "^5 + " + Str(x3) + "^5 = " + Str(y) + "^5") Goto finish EndIf Next x3 Next x2 Next x1 Next x0 PrintN("No solution was found") finish: PrintN("") PrintN("Press any key to close the console") Repeat: Delay(10) : Until Inkey() <> "" CloseConsole()
EndIf </lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
Python
Procedural
<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)
Composition of pure functions
<lang python>Euler's sum of powers conjecture
from itertools import (chain, takewhile)
- main :: IO ()
def main():
Search for counter-example
xs = enumFromTo(1)(249)
powerMap = {x**5: x for x in xs} sumMap = { x**5 + y**5: (x, y) for x in xs[1:] for y in xs if x > y }
# isExample :: (Int, Int) -> Bool def isExample(ps): p, s = ps return p - s in sumMap
# display :: (Int, Int) -> String def display(ps): p, s = ps a, b = sumMap[p - s] c, d = sumMap[s] return '^5 + '.join([str(n) for n in [a, b, c, d]]) + ( '^5 = ' + str(powerMap[p]) + '^5' )
print(__doc__ + ' – counter-example:\n') print( maybe('No counter-example found.')(display)( find(isExample)( bind(powerMap.keys())( lambda p: bind( takewhile( lambda x: p > x, sumMap.keys() ) )(lambda s: [(p, s)]) ) ) ) )
- ----------------------- GENERIC ------------------------
- Just :: a -> Maybe a
def Just(x):
Constructor for an inhabited Maybe (option type) value. Wrapper containing the result of a computation. return {'type': 'Maybe', 'Nothing': False, 'Just': x}
- Nothing :: () -> Maybe a
def Nothing():
Constructor for an empty Maybe (option type) value. Empty wrapper returned where a computation is not possible. return {'type': 'Maybe', 'Nothing': True}
- bind (>>=) :: [a] -> (a -> [b]) -> [b]
def bind(xs):
List monad injection operator. Two computations sequentially composed, with any value produced by the first passed as an argument to the second. def go(f): return chain.from_iterable(map(f, xs)) return go
- enumFromTo :: (Int, Int) -> [Int]
def enumFromTo(m):
Integer enumeration from m to n. return lambda n: range(m, 1 + n)
- find :: (a -> Bool) -> [a] -> Maybe a
def find(p):
Just the first element in the list that matches p, or Nothing if no elements match. def go(xs): try: return Just(next(x for x in xs if p(x))) except StopIteration: return Nothing() return go
- maybe :: b -> (a -> b) -> Maybe a -> b
def maybe(v):
Either the default value v, if m is Nothing, or the application of f to x, where m is Just(x). return lambda f: lambda m: v if ( None is m or m.get('Nothing') ) else f(m.get('Just'))
- MAIN ---
if __name__ == '__main__':
main()</lang>
- Output:
Euler's sum of powers conjecture – counter-example: 133^5 + 110^5 + 84^5 + 27^5 = 144^5
QL SuperBASIC
This program enhances a modular brute-force search, posted on fidonet in the 1980s, via number theoretic enhancements as used by the program listed under ZX Spectrum Basic, but without the early exit in the control framework so as to be backward-compatible with the lack of floating point in Sinclair ZX80 BASIC (whereby the latter is truly 'zeroeth' generation). To emulate running on a ZX80 (needing a 16K RAM pack & MODifications, some given below) and complete the task sooner than that for the OEM Spectrum, it relies entirely on integer functions, their upper limit being 2^15- 1 in ZX80 BASIC as well. Thus, the "slide rule" calculation of each percentage on the Spectrum is replaced by that of ones' digits across "abaci" of relatively prime bases Pi. Given that each Pi is to be <= 2^7 for said limit's sake, it takes six prime numbers or powers thereof to serve as bases of such a mixed-base number system, since it is necessary that ΠPi > 249^5 for unambiguous representation (as character strings). On ZX80s there are at most 64 consecutive printable characters (in the inverse video block: t%=48 thus becomes T=128). Just seven bases Pi <= 2^6 will be needed when the difference between 64 & a base is expressible as a four-bit offset, by which one must 'multiply' (since Z80s lack MUL) in the reduction step of the optimal assembly algorithm for MOD Pi. Such bases are: 49, 53, 55, 57, 59, 61, 64. In disproving Euler's conjecture, the program demonstrates that using 60 bits of integer precision in 1966 was 2-fold overkill, or even more so in terms of overhead cost vis-a-vis contemporaneous computers less sophisticated than a CDC 6600.
<lang qbasic> 1 CLS 2 DIM i%(255,6) : DIM a%(6) : DIM u%(6) 3 DIM v%(255,6) : DIM b%(6) : DIM d%(29) 4 RESTORE 137 6 FOR m=0 TO 6 7 READ t% 8 FOR j=1 TO 255 11 LET i%(j,m)=j MOD t% 12 LET v%(j,m)=(i%(j,m) * i%(j,m))MOD t% 14 LET v%(j,m)=(v%(j,m) * v%(j,m))MOD t% 15 LET v%(j,m)=(v%(j,m) * i%(j,m))MOD t% 17 END FOR j : END FOR m 19 PRINT "Abaci ready" 21 FOR j=10 TO 29: d%(j)=210+ j 24 FOR j=0 TO 9: d%(j)=240+ j 25 LET t%=48 30 FOR w=6 TO 246 STEP 3 33 LET o=w 42 FOR x=4 TO 248 STEP 2 44 IF o<x THEN o=x 46 FOR m=1 TO 6: a%(m)=i%((v%(w,m)+v%(x,m)),m) 50 FOR y=10 TO 245 STEP 5 54 IF o<y THEN o=y 56 FOR m=1 TO 6: b%(m)=i%((a%(m)+v%(y,m)),m) 57 FOR z=14 TO 245 STEP 7 59 IF o<z THEN o=z 60 FOR m=1 TO 6: u%(m)=i%((b%(m)+v%(z,m)),m) 65 LET s$="" : FOR m=1 TO 6: s$=s$&CHR$(u%(m)+t%) 70 LET o=o+1 : j=d%(i%((i%(w,0)+i%(x,0)+i%(y,0)+i%(z,0)),0)) 80 FOR k=j TO o STEP -30 85 LET e$="" : FOR m=1 TO 6: e$=e$&CHR$(v%(k,m)+t%) 90 IF s$=e$ THEN PRINT w,x,y,z,k,s$,e$: STOP 95 END FOR k : END FOR z : END FOR y : END FOR x : END FOR w 137 DATA 30,97,113,121,125,127,128 </lang>
- Output:
Abaci ready
27 84 110 133 144 bT`íα0 bT`íα0
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 (inexact->exact (round (expt sum 1/5))))) (break))))</lang>
- Output:
(133 110 84 27 144)
Raku
(formerly Perl 6)
<lang perl6>constant MAX = 250;
my %po5{Int}; my %sum2{Int};
(1..MAX).map: -> $i {
%po5{$i**5} = $i; for 1..MAX -> $j { %sum2{$i**5 + $j**5} = ($i, $j); }
}
my @sk = %sum2.keys.sort; %po5.keys.sort.race.map: -> $p {
for @sk -> $s { next if $p <= $s; if %sum2{$p - $s} { say ((sort |%sum2{$s}[],|%sum2{$p-$s}[]) X~ '⁵').join(' + ') ~ " = %po5{$p}" ~ "⁵"; exit; } }
}</lang>
- Output:
27⁵ + 84⁵ + 110⁵ + 133⁵ = 144⁵
REXX
fast computation
Programming note: the 3rd argument can be specified which causes an attempt to find N solutions.
The starting and ending (low and high) values can also be specified (to limit or expand the search range).
If any of the arguments are omitted, they default to the Rosetta Code task's specifications.
The method used is:
- precompute all powers of five (within the confines of allowed integers)
- precompute all (positive) differences between two applicable 5th powers
- 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 user ID G. Brougnard}
- see if the sum of any four 5th powers is equal to any 5th power
- (this is needed as the fourth number d isn't known yet).
- {all of the above utilizes REXX's sparse stemmed array hashing which eliminates the need for sorting.}
By implementing (user ID) G. Brougnard's idea of differences of two 5th powers,
the time used for computation was reduced by over a factor of seventy.
In essence, the new formula being solved is: a5 + b5 + c5 == x5 ─ d5
which lends itself to algorithm optimization by (only) having to:
- [the right side of the above equation] pre-compute all possible differences between any two applicable
integer powers of five (there are 30,135 unique differences) - [the left side of the above equation] sum any applicable three integer powers of five
- [the == part of the above equation] see if any of the above left─side sums match any of the ≈30k right─side differences
- [the right side of the above equation] pre-compute all possible differences between any two applicable
<lang rexx>/*REXX program finds unique positive integers for ────────── 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 /* " " " " " " */ w= length(H) /*W: used for display aligned numbers.*/ say center(' 'subword(sourceLine(1), 9, 3)" ", 70 +5*w, '─') /*show title from 1st line*/ numeric digits 1000 /*be able to handle the next expression*/ numeric digits max(9, length(3*H**5) ) /* " " " " 3* [H to 5th power]*/ bH= H - 2; cH= H - 1 /*calculate the upper DO loop limits.*/ !.= 0 /* [↓] define values of 5th powers. */
do pow=1 for H; @.pow= pow**5; _= @.pow; !._= 1; $._= pow end /*pow*/
?.= !.
do j=4 for H-3 /*use the range of: four to cH. */ do k=j+1 to H; _= @.k - @.j; ?._= 1 /*compute the xⁿ - dⁿ differences.*/ end /*k*/ /* [↑] diff. is always positive as k>j*/ end /*j*/ /*define [↑] 5th power differences.*/
- = 0 /*#: is the number of solutions found.*/ /* [↓] for N=∞ solutions.*/
do a=L to H-3 /*traipse through possible A values. */ /*◄──done 246 times.*/ do b=a+1 to bH; s1= @.a + @.b /* " " " B " */ /*◄──done 30,381 times.*/ do c=b+1 to cH; s2= s1 + @.c /* " " " C " */ /*◄──done 2,511,496 times.*/ if ?.s2 then do d=c+1 to H; s= s2+@.d /*find the appropriate solution. */ if !.s then call show /*Is it a solution? Then display it. */ end /*d*/ /* [↑] !.S is a boolean. */ end /*c*/ end /*b*/ end /*a*/
if #==0 then say "Didn't find a solution."; exit 0 /*──────────────────────────────────────────────────────────────────────────────────────*/ show: _= left(, 5); #= # + 1 /*_: used as a spacer; bump # counter.*/
say _ 'solution' right(#, length(N))":" _ 'a='right(a, w) _ "b="right(b, w), _ 'c='right(c, w) _ "d="right(d, w) _ 'x='right($.s, w+1) if #<N then return /*return, keep searching for more sols.*/ exit # /*stick a fork in it, we're all done. */
</lang>
- output when using the default input:
──────────────────────────── aⁿ+bⁿ+cⁿ+dⁿ==xⁿ where ⁿ=5 ───────────────────────────── solution 1: a= 27 b= 84 c=110 d=133 x= 144
- output when using the input of: 1 4000 999
─────────────────────────────── aⁿ+bⁿ+cⁿ+dⁿ==xⁿ where ⁿ=5 ─────────────────────────────── solution 1: a= 27 b= 84 c= 110 d= 133 x= 144 solution 2: a= 54 b= 168 c= 220 d= 266 x= 288 solution 3: a= 81 b= 252 c= 330 d= 399 x= 432 solution 4: a= 108 b= 336 c= 440 d= 532 x= 576 solution 5: a= 135 b= 420 c= 550 d= 665 x= 720 solution 6: a= 162 b= 504 c= 660 d= 798 x= 864 solution 7: a= 189 b= 588 c= 770 d= 931 x= 1008 solution 8: a= 216 b= 672 c= 880 d=1064 x= 1152 solution 9: a= 243 b= 756 c= 990 d=1197 x= 1296 solution 10: a= 270 b= 840 c=1100 d=1330 x= 1440 solution 11: a= 297 b= 924 c=1210 d=1463 x= 1584 solution 12: a= 324 b=1008 c=1320 d=1596 x= 1728 solution 13: a= 351 b=1092 c=1430 d=1729 x= 1872 solution 14: a= 378 b=1176 c=1540 d=1862 x= 2016 solution 15: a= 405 b=1260 c=1650 d=1995 x= 2160 solution 16: a= 432 b=1344 c=1760 d=2128 x= 2304 solution 17: a= 459 b=1428 c=1870 d=2261 x= 2448 solution 18: a= 486 b=1512 c=1980 d=2394 x= 2592 solution 19: a= 513 b=1596 c=2090 d=2527 x= 2736 solution 20: a= 540 b=1680 c=2200 d=2660 x= 2880 solution 21: a= 567 b=1764 c=2310 d=2793 x= 3024 solution 22: a= 594 b=1848 c=2420 d=2926 x= 3168 solution 23: a= 621 b=1932 c=2530 d=3059 x= 3312 solution 24: a= 648 b=2016 c=2640 d=3192 x= 3456 solution 25: a= 675 b=2100 c=2750 d=3325 x= 3600 solution 26: a= 702 b=2184 c=2860 d=3458 x= 3744 solution 27: a= 729 b=2268 c=2970 d=3591 x= 3888
lightning-fast computation
Programming note: it can be observed from the 1st REXX program's output (2nd example) that all of the index solutions are just multiples of the 1st known set:
for A, a multiple of 27 for B, a multiple of 84 for C, a multiple of 110 for D, a multiple of 133 for X, a multiple of 144
where "a multiple" is some positive integer.
Intrepid and resourceful Rosetta Code userid Pat Garrett has found that in a research paper that Jim Frye found another solution:
- 555 + 31835 + 289695 + 852825 = 853595
The paper can be seen at: A variety of Euler's conjecture.
So this 2nd known set was added to the program below.
So, index solutions are also multiples of the 2nd known set:
for A, a multiple of 55 for B, a multiple of 3,183 for C, a multiple of 28,969 for D, a multiple of 85,282 for X, a multiple of 85,359
Execution time for computing/displaying/writing 200 numbers on a slow PC is about 1/16 second.
Execution time on a slow PC for computing (alone) for 6,000 numbers is less than one second.
Execution time on a fast PC for computing (alone) for 23,686 numbers is exactly 1.00 seconds. <lang rexx>/*REXX program shows unique positive integers for ────────── aⁿ+bⁿ+cⁿ+dⁿ==xⁿ where n=5 */ numeric digits 1000 /*ensure enough decimal digs for powers*/ parse arg N oFID . /*obtain optional arguments from the CL*/ if N== | N=="," then N= 200 /*Not specified? Then use the default.*/ if oFID==|oFID=="," then oFID= 'EULERSUM.OUT' /* " " " " " " */ tell= N>=0 /*if N is ≥ 0, show output to terminal.*/ N= abs(N) /*use the absolute value of N. */
a.1= 27 ; a.2= 55 /*the A values for the two sets. */ b.1= 84 ; b.2= 3183 /* " B " " " " " */ c.1= 110 ; c.2= 28969 /* " C " " " " " */ d.1= 133 ; d.2= 85282 /* " D " " " " " */ x.1= 144 ; x.2= 85359 /* " X " " " " " */
w= length( commas(N * x.2) ) /*W: used to align displayed numbers. */ $= center(' 'subword( sourceLine(1), 9, 3)" ", 70 +5*w, '─') /*create a title.*/ call show /*show a title (from 1st line of pgm).*/ pad= left(,5) /*used for padding (spacing) the output*/ oo= 1; tt= 1 /*a counter for the A.1 & A.2 sets.*/
- = 0 /*count of number of solutions so far. */
do j=1 until #>N /*step through the possible solutions. */ one= a.1 * oo /*calculate the 1st set's A.1 value. */ two= a.2 * tt /* " " 2nd " A.2 " */ use= min(one, two) /*pick which "set" that is to be used. */ #= # + 1 /*bump counter for number of solutions.*/ if one==use then do; mult=oo; oo= oo + 1; which= 1; end if two==use then do; mult=tt; tt= tt + 1; which= 2; end $= pad 'solution' right(#,length(N))": " 'a='right( commas(a.which * mult), w), pad 'b='right( commas(b.which * mult), w), pad 'c='right( commas(c.which * mult), w), pad 'd='right( commas(d.which * mult), w), pad 'x='right( commas(x.which * mult), w) call show /*write; maybe show output to terminal.*/ res= (x.which * mult) **5 /*compute the sum of the right side. */ sum= (a.which * mult) **5 + , /* " " " " " left " */ (b.which * mult) **5 + , (c.which * mult) **5 + , (d.which * mult) **5 if sum==res then iterate /*All is kosher? Then keep truckin'. */ $= "***error*** the left side sum doesn't equal the right side result (X**5)." tell=1; call show; exit 13 /*force telling of error to terminal. */ end /*j*/
tell=1; call show $= pad ' Showed ' commas(N) " solutions, output written to file: " oFID; call show exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg _; do jc=length(_)-3 to 1 by -3; _=insert(',', _, jc); end; return _ show: if tell then say $; call lineout oFID, $; $=; return /*show and/or write it*/</lang>
- output when using the default input of: 200
(Shown at three-quarter size.)
────────────────────────────────────────────── aⁿ+bⁿ+cⁿ+dⁿ==xⁿ where n=5 ────────────────────────────────────────────── solution 1: a= 27 b= 84 c= 110 d= 133 x= 144 solution 2: a= 54 b= 168 c= 220 d= 266 x= 288 solution 3: a= 55 b= 3,183 c= 28,969 d= 85,282 x= 85,359 solution 4: a= 81 b= 252 c= 330 d= 399 x= 432 solution 5: a= 108 b= 336 c= 440 d= 532 x= 576 solution 6: a= 110 b= 6,366 c= 57,938 d= 170,564 x= 170,718 solution 7: a= 135 b= 420 c= 550 d= 665 x= 720 solution 8: a= 162 b= 504 c= 660 d= 798 x= 864 solution 9: a= 165 b= 9,549 c= 86,907 d= 255,846 x= 256,077 solution 10: a= 189 b= 588 c= 770 d= 931 x= 1,008 solution 11: a= 216 b= 672 c= 880 d= 1,064 x= 1,152 solution 12: a= 220 b= 12,732 c= 115,876 d= 341,128 x= 341,436 solution 13: a= 243 b= 756 c= 990 d= 1,197 x= 1,296 solution 14: a= 270 b= 840 c= 1,100 d= 1,330 x= 1,440 solution 15: a= 275 b= 15,915 c= 144,845 d= 426,410 x= 426,795 solution 16: a= 297 b= 924 c= 1,210 d= 1,463 x= 1,584 solution 17: a= 324 b= 1,008 c= 1,320 d= 1,596 x= 1,728 solution 18: a= 330 b= 19,098 c= 173,814 d= 511,692 x= 512,154 solution 19: a= 351 b= 1,092 c= 1,430 d= 1,729 x= 1,872 solution 20: a= 378 b= 1,176 c= 1,540 d= 1,862 x= 2,016 solution 21: a= 385 b= 22,281 c= 202,783 d= 596,974 x= 597,513 solution 22: a= 405 b= 1,260 c= 1,650 d= 1,995 x= 2,160 solution 23: a= 432 b= 1,344 c= 1,760 d= 2,128 x= 2,304 solution 24: a= 440 b= 25,464 c= 231,752 d= 682,256 x= 682,872 solution 25: a= 459 b= 1,428 c= 1,870 d= 2,261 x= 2,448 solution 26: a= 486 b= 1,512 c= 1,980 d= 2,394 x= 2,592 solution 27: a= 495 b= 28,647 c= 260,721 d= 767,538 x= 768,231 solution 28: a= 513 b= 1,596 c= 2,090 d= 2,527 x= 2,736 solution 29: a= 540 b= 1,680 c= 2,200 d= 2,660 x= 2,880 solution 30: a= 550 b= 31,830 c= 289,690 d= 852,820 x= 853,590 solution 31: a= 567 b= 1,764 c= 2,310 d= 2,793 x= 3,024 solution 32: a= 594 b= 1,848 c= 2,420 d= 2,926 x= 3,168 solution 33: a= 605 b= 35,013 c= 318,659 d= 938,102 x= 938,949 solution 34: a= 621 b= 1,932 c= 2,530 d= 3,059 x= 3,312 solution 35: a= 648 b= 2,016 c= 2,640 d= 3,192 x= 3,456 solution 36: a= 660 b= 38,196 c= 347,628 d= 1,023,384 x= 1,024,308 solution 37: a= 675 b= 2,100 c= 2,750 d= 3,325 x= 3,600 solution 38: a= 702 b= 2,184 c= 2,860 d= 3,458 x= 3,744 solution 39: a= 715 b= 41,379 c= 376,597 d= 1,108,666 x= 1,109,667 solution 40: a= 729 b= 2,268 c= 2,970 d= 3,591 x= 3,888 solution 41: a= 756 b= 2,352 c= 3,080 d= 3,724 x= 4,032 solution 42: a= 770 b= 44,562 c= 405,566 d= 1,193,948 x= 1,195,026 solution 43: a= 783 b= 2,436 c= 3,190 d= 3,857 x= 4,176 solution 44: a= 810 b= 2,520 c= 3,300 d= 3,990 x= 4,320 solution 45: a= 825 b= 47,745 c= 434,535 d= 1,279,230 x= 1,280,385 solution 46: a= 837 b= 2,604 c= 3,410 d= 4,123 x= 4,464 solution 47: a= 864 b= 2,688 c= 3,520 d= 4,256 x= 4,608 solution 48: a= 880 b= 50,928 c= 463,504 d= 1,364,512 x= 1,365,744 solution 49: a= 891 b= 2,772 c= 3,630 d= 4,389 x= 4,752 solution 50: a= 918 b= 2,856 c= 3,740 d= 4,522 x= 4,896 solution 51: a= 935 b= 54,111 c= 492,473 d= 1,449,794 x= 1,451,103 solution 52: a= 945 b= 2,940 c= 3,850 d= 4,655 x= 5,040 solution 53: a= 972 b= 3,024 c= 3,960 d= 4,788 x= 5,184 solution 54: a= 990 b= 57,294 c= 521,442 d= 1,535,076 x= 1,536,462 solution 55: a= 999 b= 3,108 c= 4,070 d= 4,921 x= 5,328 solution 56: a= 1,026 b= 3,192 c= 4,180 d= 5,054 x= 5,472 solution 57: a= 1,045 b= 60,477 c= 550,411 d= 1,620,358 x= 1,621,821 solution 58: a= 1,053 b= 3,276 c= 4,290 d= 5,187 x= 5,616 solution 59: a= 1,080 b= 3,360 c= 4,400 d= 5,320 x= 5,760 solution 60: a= 1,100 b= 63,660 c= 579,380 d= 1,705,640 x= 1,707,180 solution 61: a= 1,107 b= 3,444 c= 4,510 d= 5,453 x= 5,904 solution 62: a= 1,134 b= 3,528 c= 4,620 d= 5,586 x= 6,048 solution 63: a= 1,155 b= 66,843 c= 608,349 d= 1,790,922 x= 1,792,539 solution 64: a= 1,161 b= 3,612 c= 4,730 d= 5,719 x= 6,192 solution 65: a= 1,188 b= 3,696 c= 4,840 d= 5,852 x= 6,336 solution 66: a= 1,210 b= 70,026 c= 637,318 d= 1,876,204 x= 1,877,898 solution 67: a= 1,215 b= 3,780 c= 4,950 d= 5,985 x= 6,480 solution 68: a= 1,242 b= 3,864 c= 5,060 d= 6,118 x= 6,624 solution 69: a= 1,265 b= 73,209 c= 666,287 d= 1,961,486 x= 1,963,257 solution 70: a= 1,269 b= 3,948 c= 5,170 d= 6,251 x= 6,768 solution 71: a= 1,296 b= 4,032 c= 5,280 d= 6,384 x= 6,912 solution 72: a= 1,320 b= 76,392 c= 695,256 d= 2,046,768 x= 2,048,616 solution 73: a= 1,323 b= 4,116 c= 5,390 d= 6,517 x= 7,056 solution 74: a= 1,350 b= 4,200 c= 5,500 d= 6,650 x= 7,200 solution 75: a= 1,375 b= 79,575 c= 724,225 d= 2,132,050 x= 2,133,975 solution 76: a= 1,377 b= 4,284 c= 5,610 d= 6,783 x= 7,344 solution 77: a= 1,404 b= 4,368 c= 5,720 d= 6,916 x= 7,488 solution 78: a= 1,430 b= 82,758 c= 753,194 d= 2,217,332 x= 2,219,334 solution 79: a= 1,431 b= 4,452 c= 5,830 d= 7,049 x= 7,632 solution 80: a= 1,458 b= 4,536 c= 5,940 d= 7,182 x= 7,776 solution 81: a= 1,485 b= 85,941 c= 782,163 d= 2,302,614 x= 2,304,693 solution 82: a= 1,512 b= 4,704 c= 6,160 d= 7,448 x= 8,064 solution 83: a= 1,539 b= 4,788 c= 6,270 d= 7,581 x= 8,208 solution 84: a= 1,540 b= 89,124 c= 811,132 d= 2,387,896 x= 2,390,052 solution 85: a= 1,566 b= 4,872 c= 6,380 d= 7,714 x= 8,352 solution 86: a= 1,593 b= 4,956 c= 6,490 d= 7,847 x= 8,496 solution 87: a= 1,595 b= 92,307 c= 840,101 d= 2,473,178 x= 2,475,411 solution 88: a= 1,620 b= 5,040 c= 6,600 d= 7,980 x= 8,640 solution 89: a= 1,647 b= 5,124 c= 6,710 d= 8,113 x= 8,784 solution 90: a= 1,650 b= 95,490 c= 869,070 d= 2,558,460 x= 2,560,770 solution 91: a= 1,674 b= 5,208 c= 6,820 d= 8,246 x= 8,928 solution 92: a= 1,701 b= 5,292 c= 6,930 d= 8,379 x= 9,072 solution 93: a= 1,705 b= 98,673 c= 898,039 d= 2,643,742 x= 2,646,129 solution 94: a= 1,728 b= 5,376 c= 7,040 d= 8,512 x= 9,216 solution 95: a= 1,755 b= 5,460 c= 7,150 d= 8,645 x= 9,360 solution 96: a= 1,760 b= 101,856 c= 927,008 d= 2,729,024 x= 2,731,488 solution 97: a= 1,782 b= 5,544 c= 7,260 d= 8,778 x= 9,504 solution 98: a= 1,809 b= 5,628 c= 7,370 d= 8,911 x= 9,648 solution 99: a= 1,815 b= 105,039 c= 955,977 d= 2,814,306 x= 2,816,847 solution 100: a= 1,836 b= 5,712 c= 7,480 d= 9,044 x= 9,792 solution 101: a= 1,863 b= 5,796 c= 7,590 d= 9,177 x= 9,936 solution 102: a= 1,870 b= 108,222 c= 984,946 d= 2,899,588 x= 2,902,206 solution 103: a= 1,890 b= 5,880 c= 7,700 d= 9,310 x= 10,080 solution 104: a= 1,917 b= 5,964 c= 7,810 d= 9,443 x= 10,224 solution 105: a= 1,925 b= 111,405 c= 1,013,915 d= 2,984,870 x= 2,987,565 solution 106: a= 1,944 b= 6,048 c= 7,920 d= 9,576 x= 10,368 solution 107: a= 1,971 b= 6,132 c= 8,030 d= 9,709 x= 10,512 solution 108: a= 1,980 b= 114,588 c= 1,042,884 d= 3,070,152 x= 3,072,924 solution 109: a= 1,998 b= 6,216 c= 8,140 d= 9,842 x= 10,656 solution 110: a= 2,025 b= 6,300 c= 8,250 d= 9,975 x= 10,800 solution 111: a= 2,035 b= 117,771 c= 1,071,853 d= 3,155,434 x= 3,158,283 solution 112: a= 2,052 b= 6,384 c= 8,360 d= 10,108 x= 10,944 solution 113: a= 2,079 b= 6,468 c= 8,470 d= 10,241 x= 11,088 solution 114: a= 2,090 b= 120,954 c= 1,100,822 d= 3,240,716 x= 3,243,642 solution 115: a= 2,106 b= 6,552 c= 8,580 d= 10,374 x= 11,232 solution 116: a= 2,133 b= 6,636 c= 8,690 d= 10,507 x= 11,376 solution 117: a= 2,145 b= 124,137 c= 1,129,791 d= 3,325,998 x= 3,329,001 solution 118: a= 2,160 b= 6,720 c= 8,800 d= 10,640 x= 11,520 solution 119: a= 2,187 b= 6,804 c= 8,910 d= 10,773 x= 11,664 solution 120: a= 2,200 b= 127,320 c= 1,158,760 d= 3,411,280 x= 3,414,360 solution 121: a= 2,214 b= 6,888 c= 9,020 d= 10,906 x= 11,808 solution 122: a= 2,241 b= 6,972 c= 9,130 d= 11,039 x= 11,952 solution 123: a= 2,255 b= 130,503 c= 1,187,729 d= 3,496,562 x= 3,499,719 solution 124: a= 2,268 b= 7,056 c= 9,240 d= 11,172 x= 12,096 solution 125: a= 2,295 b= 7,140 c= 9,350 d= 11,305 x= 12,240 solution 126: a= 2,310 b= 133,686 c= 1,216,698 d= 3,581,844 x= 3,585,078 solution 127: a= 2,322 b= 7,224 c= 9,460 d= 11,438 x= 12,384 solution 128: a= 2,349 b= 7,308 c= 9,570 d= 11,571 x= 12,528 solution 129: a= 2,365 b= 136,869 c= 1,245,667 d= 3,667,126 x= 3,670,437 solution 130: a= 2,376 b= 7,392 c= 9,680 d= 11,704 x= 12,672 solution 131: a= 2,403 b= 7,476 c= 9,790 d= 11,837 x= 12,816 solution 132: a= 2,420 b= 140,052 c= 1,274,636 d= 3,752,408 x= 3,755,796 solution 133: a= 2,430 b= 7,560 c= 9,900 d= 11,970 x= 12,960 solution 134: a= 2,457 b= 7,644 c= 10,010 d= 12,103 x= 13,104 solution 135: a= 2,475 b= 143,235 c= 1,303,605 d= 3,837,690 x= 3,841,155 solution 136: a= 2,484 b= 7,728 c= 10,120 d= 12,236 x= 13,248 solution 137: a= 2,511 b= 7,812 c= 10,230 d= 12,369 x= 13,392 solution 138: a= 2,530 b= 146,418 c= 1,332,574 d= 3,922,972 x= 3,926,514 solution 139: a= 2,538 b= 7,896 c= 10,340 d= 12,502 x= 13,536 solution 140: a= 2,565 b= 7,980 c= 10,450 d= 12,635 x= 13,680 solution 141: a= 2,585 b= 149,601 c= 1,361,543 d= 4,008,254 x= 4,011,873 solution 142: a= 2,592 b= 8,064 c= 10,560 d= 12,768 x= 13,824 solution 143: a= 2,619 b= 8,148 c= 10,670 d= 12,901 x= 13,968 solution 144: a= 2,640 b= 152,784 c= 1,390,512 d= 4,093,536 x= 4,097,232 solution 145: a= 2,646 b= 8,232 c= 10,780 d= 13,034 x= 14,112 solution 146: a= 2,673 b= 8,316 c= 10,890 d= 13,167 x= 14,256 solution 147: a= 2,695 b= 155,967 c= 1,419,481 d= 4,178,818 x= 4,182,591 solution 148: a= 2,700 b= 8,400 c= 11,000 d= 13,300 x= 14,400 solution 149: a= 2,727 b= 8,484 c= 11,110 d= 13,433 x= 14,544 solution 150: a= 2,750 b= 159,150 c= 1,448,450 d= 4,264,100 x= 4,267,950 solution 151: a= 2,754 b= 8,568 c= 11,220 d= 13,566 x= 14,688 solution 152: a= 2,781 b= 8,652 c= 11,330 d= 13,699 x= 14,832 solution 153: a= 2,805 b= 162,333 c= 1,477,419 d= 4,349,382 x= 4,353,309 solution 154: a= 2,808 b= 8,736 c= 11,440 d= 13,832 x= 14,976 solution 155: a= 2,835 b= 8,820 c= 11,550 d= 13,965 x= 15,120 solution 156: a= 2,860 b= 165,516 c= 1,506,388 d= 4,434,664 x= 4,438,668 solution 157: a= 2,862 b= 8,904 c= 11,660 d= 14,098 x= 15,264 solution 158: a= 2,889 b= 8,988 c= 11,770 d= 14,231 x= 15,408 solution 159: a= 2,915 b= 168,699 c= 1,535,357 d= 4,519,946 x= 4,524,027 solution 160: a= 2,916 b= 9,072 c= 11,880 d= 14,364 x= 15,552 solution 161: a= 2,943 b= 9,156 c= 11,990 d= 14,497 x= 15,696 solution 162: a= 2,970 b= 171,882 c= 1,564,326 d= 4,605,228 x= 4,609,386 solution 163: a= 2,997 b= 9,324 c= 12,210 d= 14,763 x= 15,984 solution 164: a= 3,024 b= 9,408 c= 12,320 d= 14,896 x= 16,128 solution 165: a= 3,025 b= 175,065 c= 1,593,295 d= 4,690,510 x= 4,694,745 solution 166: a= 3,051 b= 9,492 c= 12,430 d= 15,029 x= 16,272 solution 167: a= 3,078 b= 9,576 c= 12,540 d= 15,162 x= 16,416 solution 168: a= 3,080 b= 178,248 c= 1,622,264 d= 4,775,792 x= 4,780,104 solution 169: a= 3,105 b= 9,660 c= 12,650 d= 15,295 x= 16,560 solution 170: a= 3,132 b= 9,744 c= 12,760 d= 15,428 x= 16,704 solution 171: a= 3,135 b= 181,431 c= 1,651,233 d= 4,861,074 x= 4,865,463 solution 172: a= 3,159 b= 9,828 c= 12,870 d= 15,561 x= 16,848 solution 173: a= 3,186 b= 9,912 c= 12,980 d= 15,694 x= 16,992 solution 174: a= 3,190 b= 184,614 c= 1,680,202 d= 4,946,356 x= 4,950,822 solution 175: a= 3,213 b= 9,996 c= 13,090 d= 15,827 x= 17,136 solution 176: a= 3,240 b= 10,080 c= 13,200 d= 15,960 x= 17,280 solution 177: a= 3,245 b= 187,797 c= 1,709,171 d= 5,031,638 x= 5,036,181 solution 178: a= 3,267 b= 10,164 c= 13,310 d= 16,093 x= 17,424 solution 179: a= 3,294 b= 10,248 c= 13,420 d= 16,226 x= 17,568 solution 180: a= 3,300 b= 190,980 c= 1,738,140 d= 5,116,920 x= 5,121,540 solution 181: a= 3,321 b= 10,332 c= 13,530 d= 16,359 x= 17,712 solution 182: a= 3,348 b= 10,416 c= 13,640 d= 16,492 x= 17,856 solution 183: a= 3,355 b= 194,163 c= 1,767,109 d= 5,202,202 x= 5,206,899 solution 184: a= 3,375 b= 10,500 c= 13,750 d= 16,625 x= 18,000 solution 185: a= 3,402 b= 10,584 c= 13,860 d= 16,758 x= 18,144 solution 186: a= 3,410 b= 197,346 c= 1,796,078 d= 5,287,484 x= 5,292,258 solution 187: a= 3,429 b= 10,668 c= 13,970 d= 16,891 x= 18,288 solution 188: a= 3,456 b= 10,752 c= 14,080 d= 17,024 x= 18,432 solution 189: a= 3,465 b= 200,529 c= 1,825,047 d= 5,372,766 x= 5,377,617 solution 190: a= 3,483 b= 10,836 c= 14,190 d= 17,157 x= 18,576 solution 191: a= 3,510 b= 10,920 c= 14,300 d= 17,290 x= 18,720 solution 192: a= 3,520 b= 203,712 c= 1,854,016 d= 5,458,048 x= 5,462,976 solution 193: a= 3,537 b= 11,004 c= 14,410 d= 17,423 x= 18,864 solution 194: a= 3,564 b= 11,088 c= 14,520 d= 17,556 x= 19,008 solution 195: a= 3,575 b= 206,895 c= 1,882,985 d= 5,543,330 x= 5,548,335 solution 196: a= 3,591 b= 11,172 c= 14,630 d= 17,689 x= 19,152 solution 197: a= 3,618 b= 11,256 c= 14,740 d= 17,822 x= 19,296 solution 198: a= 3,630 b= 210,078 c= 1,911,954 d= 5,628,612 x= 5,633,694 solution 199: a= 3,645 b= 11,340 c= 14,850 d= 17,955 x= 19,440 solution 200: a= 3,672 b= 11,424 c= 14,960 d= 18,088 x= 19,584 Showed 200 solutions, output written to file: EULERSUM.OUT
Ring
<lang ring>
- Project : Euler's sum of powers conjecture
max=250 for w = 1 to max
for x = 1 to w for y = 1 to x for z = 1 to y sum = pow(w,5) + pow(x,5) + pow(y,5) + pow(z,5) s1 = floor(pow(sum,0.2)) if sum = pow(s1,5) see "" + w + "^5 + " + x + "^5 + " + y + "^5 + " + z + "^5 = " + s1 + "^5" ok next next next
next </lang> Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
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.
Run BASIC
<lang runbasic> max=250 FOR w = 1 TO max
FOR x = 1 TO w FOR y = 1 TO x FOR z = 1 TO y sum = w^5 + x^5 + y^5 + z^5 s1 = INT(sum^0.2) IF sum=s1^5 THEN PRINT w;"^5 + ";x;"^5 + ";y;"^5 + ";z;"^5 = ";s1;"^5" end end if NEXT z NEXT y NEXT x
NEXT w</lang>
133^5 + 110^5 + 84^5 + 27^5 = 144^5
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
Scala
Functional programming
<lang Scala>import scala.collection.Searching.{Found, search}
object EulerSopConjecture extends App {
val (maxNumber, fifth) = (250, (1 to 250).map { i => math.pow(i, 5).toLong })
def binSearch(fact: Int*) = fifth.search(fact.map(f => fifth(f)).sum)
def sop = (0 until maxNumber) .flatMap(a => (a until maxNumber) .flatMap(b => (b until maxNumber) .flatMap(c => (c until maxNumber) .map { case x$1@d => (binSearch(a, b, c, d), x$1) } .withFilter { case (f, _) => f.isInstanceOf[Found] } .map { case (f, d) => (a + 1, b + 1, c + 1, d + 1, f.insertionPoint + 1) }))).take(1) .map { case (a, b, c, d, f) => s"$a⁵ + $b⁵ + $c⁵ + $d⁵ = $f⁵" }
println(sop)
}</lang>
- Output:
Vector(27⁵ + 84⁵ + 110⁵ + 133⁵ = 144⁵)
Seed7
<lang seed7>$ include "seed7_05.s7i";
const func integer: binarySearch (in array integer: arr, in integer: aKey) is func
result var integer: index is 0; local var integer: low is 1; var integer: high is 0; var integer: middle is 0; begin high := length(arr); while index = 0 and low <= high do middle := (low + high) div 2; if aKey < arr[middle] then high := pred(middle); elsif aKey > arr[middle] then low := succ(middle); else index := middle; end if; end while; end func;
const proc: main is func
local var array integer: p5 is 249 times 0; var integer: i is 0; var integer: x0 is 0; var integer: x1 is 0; var integer: x2 is 0; var integer: x3 is 0; var integer: sum is 0; var integer: y is 0; var boolean: found is FALSE; begin for i range 1 to 249 do p5[i] := i ** 5; end for; for x0 range 1 to 249 until found do for x1 range 1 to pred(x0) until found do for x2 range 1 to pred(x1) until found do for x3 range 1 to pred(x2) until found do sum := p5[x0] + p5[x1] + p5[x2] + p5[x3]; y := binarySearch(p5, sum); if y > 0 then writeln(x0 <& "**5 + " <& x1 <& "**5 + " <& x2 <& "**5 + " <& x3 <& "**5 = " <& y <& "**5"); found := TRUE; end if; end for; end for; end for; end for; if not found then writeln("No solution was found"); end if; end func;</lang>
- Output:
133**5 + 110**5 + 84**5 + 27**5 = 144**5
SenseTalk
<lang sensetalk>findEulerSumOfPowers to findEulerSumOfPowers
set MAX_NUMBER to 250 set possibleValues to 1..MAX_NUMBER set possible5thPowers to each item of possibleValues to the power of 5 repeat for x0 in 1..250 repeat for x1 in 1..x0 repeat for x2 in 1..x1 repeat for x3 in 1..x2 set possibleSum to item x0 of possible5thPowers \ plus item x1 of possible5thPowers \ plus item x2 of possible5thPowers \ plus item x3 of possible5thPowers if possibleSum is in possible5thPowers put x0 & "^5 + " & x1 & "^5 + " & x2 & "^5 + " & x3 & "^5 = " & the item number of possibleSum within possible5thPowers & "^5" return end if end repeat end repeat end repeat end repeat
end findEulerSumOfPowers</lang>
Sidef
<lang ruby>define range = (1 ..^ 250)
var p5 = Hash() var sum2 = Hash()
for i in (range) {
p5{i**5} = i for j in (range) { sum2{i**5 + j**5} = [i, j] }
}
var sk = sum2.keys.map{ Num(_) }.sort
for p in (p5.keys.map{ Num(_) }.sort) {
var s = sk.first {|s| p > s && sum2.exists(p-s) } \\ next
var t = (sum2{s} + sum2{p-s} -> map{|n| "#{n}⁵" }.join(' + ')) say "#{t} = #{p5{p}}⁵" break
}</lang>
- Output:
84⁵ + 27⁵ + 133⁵ + 110⁵ = 144⁵
Swift
<lang swift>extension BinaryInteger {
@inlinable public func power(_ n: Self) -> Self { return stride(from: 0, to: n, by: 1).lazy.map({_ in self }).reduce(1, *) }
}
func sumOfPowers(maxN: Int = 250) -> (Int, Int, Int, Int, Int) {
let pow5 = (0..<maxN).map({ $0.power(5) }) let pow5ToN = {n in pow5.firstIndex(of: n)}
for x0 in 1..<maxN { for x1 in 1..<x0 { for x2 in 1..<x1 { for x3 in 1..<x2 { let powSum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3]
if let idx = pow5ToN(powSum) { return (x0, x1, x2, x3, idx) } } } } }
fatalError("Did not find solution")
}
let (x0, x1, x2, x3, y) = sumOfPowers()
print("\(x0)^5 + \(x1)^5 + \(x2)^5 \(x3)^5 = \(y)^5")</lang>
- Output:
133^5 + 110^5 + 84^5 + 27^5 = 144^5
Tcl
<lang Tcl>proc doit {{badx 250} {complete 0}} {
## NB: $badx is exclusive upper limit, and also limits y! for {set y 1} {$y < $badx} {incr y} { set s [expr {$y ** 5}] set r5($s) $y ;# fifth roots of valid sums } for {set a 1} {$a < $badx} {incr a} { set suma [expr {$a ** 5}] for {set b 1} {$b <= $a} {incr b} { set sumb [expr {$suma + ($b ** 5)}] for {set c 1} {$c <= $b} {incr c} { set sumc [expr {$sumb + ($c ** 5)}] for {set d 1} {$d <= $c} {incr d} { set sumd [expr {$sumc + ($d ** 5)}] if {[info exists r5($sumd)]} { set e $r5($sumd) puts "$e^5 = $a^5 + $b^5 + $c^5 + $d^5" if {!$complete} { return } } } } } } puts "search complete (x < $badx)"
} doit</lang>
- Output:
144^5 = 133^5 + 110^5 + 84^5 + 27^5
real 0m2.387s
UNIX Shell
Shell is not the go-to language for number-crunching, but if you're going to use a shell, it looks like ksh is the fastest option, at about 8x faster than bash and 2x faster than zsh.
<lang bash>MAX=250 pow5=() for (( i=1; i<MAX; ++i )); do
pow5[i]=$(( i*i*i*i*i ))
done for (( a=1; a<MAX; ++a )); do
for (( b=a+1; b<MAX; ++b )); do for (( c=b+1; c<MAX; ++c )); do for (( d=c+1; d<MAX; ++d )); do (( sum=pow5[a]+pow5[b]+pow5[c]+pow5[d] )) (( low=d+3 )) (( high=MAX )) while (( low <= high )); do (( guess=(low+high)/2 )) if (( pow5[guess] == sum )); then printf 'Found example: %d⁵+%d⁵+%d⁵+%d⁵=%d⁵\n' "$a" "$b" "$c" "$d" "$guess" exit 0 elif (( pow5[guess] < sum )); then (( low=guess+1 )) else (( high=guess-1 )) fi done done done done
done printf 'No examples found.\n' exit 1</lang>
- Output:
$ time bash esop.sh Found example: 27⁵+84⁵+110⁵+133⁵=144⁵ bash esop.sh 6953.75s user 37.53s system 99% cpu 1:57:02.41 total $ time ksh esop.sh Found example: 27⁵+84⁵+110⁵+133⁵=144⁵ ksh esop.sh 855.66s user 5.30s system 99% cpu 14:26.78 total $ time zsh esop.sh Found example: 27⁵+84⁵+110⁵+133⁵=144⁵ zsh esop.sh 1969.48s user 250.82s system 99% cpu 37:11.62 total
VBA
<lang vb>Public Declare Function GetTickCount Lib "kernel32.dll" () As Long
Public Sub begin()
start_int = GetTickCount() main Debug.Print (GetTickCount() - start_int) / 1000 & " seconds"
End Sub Private Function pow(x, y) As Variant
pow = CDec(Application.WorksheetFunction.Power(x, y))
End Function Private Sub main()
For x0 = 1 To 250 For x1 = 1 To x0 For x2 = 1 To x1 For x3 = 1 To x2 sum = CDec(pow(x0, 5) + pow(x1, 5) + pow(x2, 5) + pow(x3, 5)) s1 = Int(pow(sum, 0.2)) If sum = pow(s1, 5) Then Debug.Print x0 & "^5 + " & x1 & "^5 + " & x2 & "^5 + " & x3 & "^5 = " & s1 Exit Sub End If Next x3 Next x2 Next x1 Next x0
End Sub</lang>
- Output:
33^5 + 110^5 + 84^5 + 27^5 = 144 160,187 seconds
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
Visual Basic .NET
Paired Powers Algorithm
<lang vbnet>Module Module1
Structure Pair Dim a, b As Integer Sub New(x as integer, y as integer) a = x : b = y End Sub End Structure
Dim max As Integer = 250 Dim p5() As Long, sum2 As SortedDictionary(Of Long, Pair) = New SortedDictionary(Of Long, Pair)
Function Fmt(p As Pair) As String Return String.Format("{0}^5 + {1}^5", p.a, p.b) End Function
Sub Init() p5(0) = 0 : p5(1) = 1 : For i As Integer = 1 To max - 1 For j As Integer = i + 1 To max p5(j) = CLng(j) * j : p5(j) *= p5(j) * j sum2.Add(p5(i) + p5(j), New Pair(i, j)) Next Next End Sub
Sub Calc(Optional findLowest As Boolean = True) For i As Integer = 1 To max : Dim p As Long = p5(i) For Each s In sum2.Keys Dim t As Long = p - s : If t <= 0 Then Exit For If sum2.Keys.Contains(t) AndAlso sum2.Item(t).a > sum2.Item(s).b Then Console.WriteLine(" {1} + {2} = {0}^5", i, Fmt(sum2.Item(s)), Fmt(sum2.Item(t))) If findLowest Then Exit Sub End If Next : Next End Sub
Sub Main(args As String()) If args.Count > 0 Then Dim t As Integer = 0 : Integer.TryParse(args(0), t) If t > 0 AndAlso t < 5405 Then max = t End If Console.WriteLine("Checking from 1 to {0}...", max) For i As Integer = 0 To 1 ReDim p5(max) : sum2.Clear() Dim st As DateTime = DateTime.Now Init() : Calc(i = 0) Console.WriteLine("{0} Computation time to {2} was {1} seconds{0}", vbLf, (DateTime.Now - st).TotalSeconds, If(i = 0, "find lowest one", "check entire space")) Next If Diagnostics.Debugger.IsAttached Then Console.ReadKey() End Sub
End Module</lang>
- Output:
(No command line arguments)
Checking from 1 to 250... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 Computation time to find lowest one was 0.0807819 seconds 27^5 + 84^5 + 110^5 + 133^5 = 144^5 Computation time to check entire space was 0.3830103 seconds
Command line argument = "1000"
Checking from 1 to 1000... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 Computation time to find lowest one was 0.3112007 seconds 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time to check entire space was 28.8847393 seconds
Paired Powers w/ Mod 30 Shortcut and Threading
If one divides the searched array of powers (sum2m()) into 30 pieces, the search time can be reduced by only searching the appropriate one (determined by the Mod 30 value of the value being sought). Once broken down by this, it is now easier to use threading to further reduce the computation time.
The following compares the plain paired powers algorithm to the plain powers plus the mod 30 shortcut algorithm, without and with threading.
<lang vbnet>Module Module1
Structure Pair Dim a, b As Integer Sub New(x As Integer, y As Integer) a = x : b = y End Sub End Structure
Dim min As Integer = 1, max As Integer = 250 Dim p5() As Long, sum2 As SortedDictionary(Of Long, Pair) = New SortedDictionary(Of Long, Pair), sum2m(29) As SortedDictionary(Of Long, Pair)
Function Fmt(p As Pair) As String Return String.Format("{0}^5 + {1}^5", p.a, p.b) End Function
Sub Init() p5(0) = 0 : p5(min) = CLng(min) * min : p5(min) *= p5(min) * min For i As Integer = min To max - 1 For j As Integer = i + 1 To max p5(j) = CLng(j) * j : p5(j) *= p5(j) * j If j = max Then Continue For sum2.Add(p5(i) + p5(j), New Pair(i, j)) Next Next End Sub
Sub InitM() For i As Integer = 0 To 29 : sum2m(i) = New SortedDictionary(Of Long, Pair) : Next p5(0) = 0 : p5(min) = CLng(min) * min : p5(min) *= p5(min) * min For i As Integer = min To max - 1 For j As Integer = i + 1 To max p5(j) = CLng(j) * j : p5(j) *= p5(j) * j If j = max Then Continue For Dim x As Long = p5(i) + p5(j) sum2m(x Mod 30).Add(x, New Pair(i, j)) Next Next End Sub
Sub Calc(Optional findLowest As Boolean = True) For i As Integer = min To max : Dim p As Long = p5(i) For Each s In sum2.Keys Dim t As Long = p - s : If t <= 0 Then Exit For If sum2.Keys.Contains(t) AndAlso sum2.Item(t).a > sum2.Item(s).b Then Console.WriteLine(" {1} + {2} = {0}^5", i, Fmt(sum2.Item(s)), Fmt(sum2.Item(t))) If findLowest Then Exit Sub End If Next : Next End Sub
Function CalcM(m As Integer) As List(Of String) Dim res As New List(Of String) For i As Integer = min To max Dim pm As Integer = i Mod 30, mp As Integer = (pm - m + 30) Mod 30 For Each s In sum2m(m).Keys Dim t As Long = p5(i) - s : If t <= 0 Then Exit For If sum2m(mp).Keys.Contains(t) AndAlso sum2m(mp).Item(t).a > sum2m(m).Item(s).b Then res.Add(String.Format(" {1} + {2} = {0}^5", i, Fmt(sum2m(m).Item(s)), Fmt(sum2m(mp).Item(t)))) End If Next : Next Return res End Function
Function Snip(s As String) As Integer Dim p As Integer = s.IndexOf("=") + 1 Return s.Substring(p, s.IndexOf("^", p) - p) End Function
Function CompareRes(ByVal x As String, ByVal y As String) As Integer CompareRes = Snip(x).CompareTo(Snip(y)) If CompareRes = 0 Then CompareRes = x.CompareTo(y) End Function
Function Validify(def As Integer, s As String) As Integer Validify = def : Dim t As Integer = 0 : Integer.TryParse(s, t) If t >= 1 AndAlso Math.Pow(t, 5) < (Long.MaxValue >> 1) Then Validify = t End Function
Sub Switch(ByRef a As Integer, ByRef b As Integer) Dim t As Integer = a : a = b : b = t End Sub
Sub Main(args As String()) Select Case args.Count Case 1 : max = Validify(max, args(0)) Case > 1 min = Validify(min, args(0)) max = Validify(max, args(1)) If max < min Then Switch(max, min) End Select Console.WriteLine("Paired powers, checking from {0} to {1}...", min, max) For i As Integer = 0 To 1 ReDim p5(max) : sum2.Clear() Dim st As DateTime = DateTime.Now Init() : Calc(i = 0) Console.WriteLine("{0} Computation time to {2} was {1} seconds{0}", vbLf, (DateTime.Now - st).TotalSeconds, If(i = 0, "find lowest one", "check entire space")) Next For i As Integer = 0 To 1 Console.WriteLine("Paired powers with Mod 30 shortcut (entire space) {2}, checking from {0} to {1}...", min, max, If(i = 0, "sequential", "parallel")) ReDim p5(max) Dim res As List(Of String) = New List(Of String) Dim st As DateTime = DateTime.Now Dim taskList As New List(Of Task(Of List(Of String))) InitM() Select Case i Case 0 For j As Integer = 0 To 29 res.AddRange(CalcM(j)) Next Case 1 For j As Integer = 0 To 29 : Dim jj = j taskList.Add(Task.Run(Function() CalcM(jj))) Next Task.WhenAll(taskList) For Each item In taskList.Select(Function(t) t.Result) res.AddRange(item) : Next End Select res.Sort(AddressOf CompareRes) For Each item In res Console.WriteLine(item) : Next Console.WriteLine("{0} Computation time was {1} seconds{0}", vbLf, (DateTime.Now - st).TotalSeconds) Next If Diagnostics.Debugger.IsAttached Then Console.ReadKey() End Sub
End Module</lang>
- Output:
(No command line arguments)
Paired powers, checking from 1 to 250...27^5 + 84^5 + 110^5 + 133^5 = 144^5
Computation time to find lowest one was 0.0781252 seconds
27^5 + 84^5 + 110^5 + 133^5 = 144^5
Computation time to check entire space was 0.3280574 seconds
Paired powers with Mod 30 shortcut (entire space) sequential, checking from 1 to 250... 27^5 + 84^5 + 110^5 + 133^5 = 144^5
Computation time was 0.2655529 seconds
Paired powers with Mod 30 shortcut (entire space) parallel, checking from 1 to 250... 27^5 + 84^5 + 110^5 + 133^5 = 144^5
Computation time was 0.0624651 seconds
(command line argument = "1000")
Paired powers, checking from 1 to 1000... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 Computation time to find lowest one was 0.2499343 seconds 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time to check entire space was 27.805961 seconds Paired powers with Mod 30 shortcut (entire space) sequential, checking from 1 to 1000... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time was 23.8068928 seconds Paired powers with Mod 30 shortcut (entire space) parallel, checking from 1 to 1000... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time was 5.4205943 seconds
(command line arguments = "27 864")
Paired powers, checking from 27 to 864... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 Computation time to find lowest one was 0.1562309 seconds 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time to check entire space was 15.8243802 seconds Paired powers with Mod 30 shortcut (entire space) sequential, checking from 27 to 864... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time was 13.0438215 seconds Paired powers with Mod 30 shortcut (entire space) parallel, checking from 27 to 864... 27^5 + 84^5 + 110^5 + 133^5 = 144^5 54^5 + 168^5 + 220^5 + 266^5 = 288^5 81^5 + 252^5 + 330^5 + 399^5 = 432^5 108^5 + 336^5 + 440^5 + 532^5 = 576^5 135^5 + 420^5 + 550^5 + 665^5 = 720^5 162^5 + 504^5 + 660^5 + 798^5 = 864^5 Computation time was 3.0305365 seconds
(command line arguments = "189 1008")
Paired powers, checking from 189 to 1008... 189^5 + 588^5 + 770^5 + 931^5 = 1008^5 Computation time to find lowest one was 14.6840411 seconds 189^5 + 588^5 + 770^5 + 931^5 = 1008^5 Computation time to check entire space was 14.7777685 seconds Paired powers with Mod 30 shortcut (entire space) sequential, checking from 189 to 1008... 189^5 + 588^5 + 770^5 + 931^5 = 1008^5 Computation time was 12.4814705 seconds Paired powers with Mod 30 shortcut (entire space) parallel, checking from 189 to 1008... 189^5 + 588^5 + 770^5 + 931^5 = 1008^5 Computation time was 2.7180777 seconds
Wren
<lang ecmascript>var start = System.clock var n = 250 var m = 30
var p5 = List.filled(n+m+1, 0) var s = 0 while (s < n) {
var sq = s * s p5[s] = sq * sq * s s = s + 1
} var max = p5[n-1] while (s < p5.count) {
p5[s] = max + 1 s = s + 1
} for (a in 1...n-3) {
for (b in a + 1...n-2) { for (c in b + 1...n-1) { var d = c + 1 var t = p5[a] + p5[b] + p5[c] var e = d + (t % m) s = t + p5[d] while (s <= max) { e = e - m while (p5[e+m] <= s) e = e + m if (p5[e] == s) { System.print("%(a)⁵ + %(b)⁵ + %(c)⁵ + %(d)⁵ = %(e)⁵") System.print("Took %(System.clock - start) seconds") return } d = d + 1 e = e + 1 s = t + p5[d] } } }
}</lang>
- Output:
Timing is for an Intel Core i7-8565U machine running Wren 0.2.0 on Ubuntu 18.04.
27⁵ + 84⁵ + 110⁵ + 133⁵ = 144⁵ Took 6.836934 seconds
Zig
<lang Zig> const std = @import("std"); const stdout = std.io.getStdOut().outStream();
pub fn main() !void {
var pow5: [250]i64 = undefined; for (pow5) |*e, i| { const n = @intCast(i64, i); e.* = n * n * n * n * n; } var x0: u16 = 4; while (x0 < pow5.len) : (x0 += 1) { var x1: u16 = 3; while (x1 < x0) : (x1 += 1) { var x2: u16 = 2; while (x2 < x1) : (x2 += 1) { var x3: u16 = 1; while (x3 < x2) : (x3 += 1) { const sum = pow5[x0] + pow5[x1] + pow5[x2] + pow5[x3]; var y: u16 = x0 + 1; while (y < pow5.len) : (y += 1) if (sum == pow5[y]) { try stdout.print("{}⁵ + {}⁵ + {}⁵ + {}⁵ = {}⁵\n", .{ x0, x1, x2, x3, y }); return; }; } } } }
try stdout.print("Sorry, no solution found.\n", .{});
} </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)); break(4); // the foreach is actually four loops
}</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:=Dictionary(),Dictionary(); 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
ZX Spectrum Basic
This "abacus revision" reverts back to an earlier one, i.e. the "slide rule" one calculating the logarithmic 'percentage' for each of 4 summands. It also calculates their ones' digits as if on a base m abacus, that complements the other base m digits embedded in their percentages via a check sum of the ones' digits when the percentages add up to 1. Because any other argument is less than m, there is sufficient additional precision so as to not find false solutions while seeking the first primitive solution. 1 is excluded a priori for (e.g.) w, since it is well-known that the only integral points on 1=m^5-x^5-y^5-z^5 are the obvious ones (that aren't distinct\all >0). Nonetheless, 1 (as the zeroeth 'prime' Po) can be the first of 3 factors (one of the others being a power of the first 4 primes) for specifying a LHS argument. Given 4 LHS arguments each raised to a 5th power as is the RHS for which m<=ΣΠPi<2^(3+5), i=0 to 4, the LHS solution (if one exists) will consist of 4 elements of a factorial domain that are each a triplet (a prime\Po * a prime^1st\2nd power * a prime^1st\4th power) multiple having a distinct pair of the first 4 primes present & absent. Such potential solutions are generated by aiming for simplicity rather than efficiency (e.g. not generating potential solutions where each multiple is even). Two theorems' consequences intended for an even earlier revision are also applied: Fermat's Little Theorem (as proven later by Euler) whereby Xi^5 mod P = Xi mod P for the first 3 primes; and the Chinese Remainder Theorem in reverse whereby m= k- i*2*3*5, such that k is chosen to be the highest allowed m congruent mod 30 to the sum of the LHS arguments of a potential solution. Although still slow, this revision performs the given task on any ZX Spectrum, despite being limited to 32-bit precision.
<lang zxbasic> 2 CLS 5 DIM k(29): DIM q(249) 10 FOR i=4 TO 249: LET q(i)=LN i : NEXT i 11 REM enhancements for the much expanded Spectrum Next: DIM p(248,249) 12 REM FOR j=4TO 248:FOR i=15TO 249:LET p(j,i)=EXP (q(j)-q(i))*5:NEXT i:NEXT j 14 PRINT "slide rule ready" 15 FOR i=0 TO 9: LET k(i)=240+ i : NEXT i 17 FOR i=10 TO 29: LET k(i)=210+ i : NEXT i 20 FOR w=6 TO 246 STEP 3 21 LET o=w 22 FOR x=4 TO 248 STEP 2 23 IF o<x THEN LET o=x 24 FOR y=10 TO 245 STEP 5 25 IF o<y THEN LET o=y 26 FOR z=14 TO 245 STEP 7 27 IF o<z THEN LET o=z 30 LET o=o+1 : LET m=k(FN f((w+x+y+z),30)) 34 IF m<o THEN GO TO 90 40 REM LET s=p(w,m)+p(x,m)+p(y,m)+p(z,m) instead of: 42 LET s=EXP((q(w)-q(m))*5) 43 LET s=EXP((q(x)-q(m))*5)+ s 45 LET s=EXP((q(y)-q(m))*5)+ s 47 LET s=EXP((q(z)-q(m))*5)+ s 50 IF s<>1 THEN GO TO 80 52 LET a=FN f(w*w,m) : LET a=FN f(a*a*w,m) 53 LET b=FN f(x*x,m) : LET b=FN f(b*b*x,m) 55 LET c=FN f(y*y,m) : LET c=FN f(c*c*y,m) 57 LET d=FN f(z*z,m) : LET d=FN f(d*d*z,m) 60 LET u=FN f((a+b+c+d),m) 65 IF u THEN GO TO 80 73 PRINT w;"^5+";x;"^5+";y;"^5+";z;"^5=";m;"^5": STOP 80 IF s<1 THEN m=m-30 : GO TO 34 90 NEXT z: NEXT y: NEXT x: NEXT w 100 DEF FN f(e,n)=e- INT(e/n)*n </lang>
- Output:
slide rule ready 27^5+84^5+110^5+133^5=144^5
- Programming Tasks
- Solutions by Programming Task
- 11l
- 360 Assembly
- 6502 Assembly
- Ada
- ALGOL 68
- ALGOL W
- Arturo
- AWK
- Bracmat
- C
- C sharp
- C++
- Clojure
- COBOL
- Common Lisp
- D
- Delphi
- EasyLang
- EchoLisp
- Elixir
- ERRE
- F Sharp
- Factor
- Forth
- Fortran
- FreeBASIC
- Go
- Groovy
- Haskell
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- Lua
- Mathematica
- Wolfram Language
- Microsoft Small Basic
- Modula-2
- Modula-2 examples needing attention
- Examples needing attention
- Nim
- Oforth
- PARI/GP
- Pascal
- Perl
- Phix
- PHP
- PicoLisp
- PowerShell
- Prolog
- PureBasic
- Python
- QL SuperBASIC
- Racket
- Raku
- REXX
- Ring
- Ruby
- Run BASIC
- Rust
- Scala
- Seed7
- SenseTalk
- Sidef
- Swift
- Tcl
- UNIX Shell
- VBA
- VBScript
- Visual Basic .NET
- Wren
- Zig
- Zkl
- ZX Spectrum Basic