Open question in 1937…short SAS program today

John D. Cook posted a story about Hardy, Ramanujan, and Euler and discusses a conjecture in number theory from 1937. Cook says,

Euler discovered 653,518,657 = 158^4 + 59^4 = 134^4 + 133^4 and that this was the smallest [integer] known to be the sum of two fourth powers in two ways. It seems odd now to think of such questions being unresolved. Today we’d ask Hardy “What do you mean 653518657 is the smallest known example? Why didn’t you write a little program to find out whether it really is the smallest?”

Cook goes on to encourage his readers to write a program that settles the conjecture, which is obviously no longer an open problem. Here’s a SAS/IML solution. First, introduce a function that I use a lot:

proc iml;
start ExpandGrid( x, y ); 
/* Return ordered pairs on a uniform grid of points */
   Nx = nrow(x);    Ny = nrow(y);
   x = repeat(x, Ny);
   y = shape( repeat(y, 1, Nx), 0, 1 );
   return ( x || y );
finish;

The ExpandGrid function generates a uniform grid of all ordered pairs of it’s arguments. To solve the “conjecture,” generate all ordered pairs that contain the integers 1–159 and check to see if the sum of fourth powers contain duplicate values. The number 159 is the greatest integer that is less than the fourth root of 653,518,657.

/* settle conjecture of Euler/Hardy: Is 653,518,657 the smallest
   integer that can be written as the sum of 4th powers in two ways? */
k = floor(653518657##0.25);    /* only need to test up to 4th root */
g = ExpandGrid(t(1:k), t(1:k));/* all ordered pairs */
g = g[ loc(g[,1] <= g[,2]), ]; /* omit symmetric pairs */
f = g[,1]##4 + g[,2]##4;       /* sum of 4th powers */
print (nrow(f)) (ncol(unique(f))); /* exactly one sum is the same */

The output proves the conjecture: there are 12,720 ordered pairs of integers to consider, and there are 12,719 unique sums of fourth powers. That means that there is exactly one sum of fourth powers that can be written in two different ways, and Euler has already told us the particular values.

If you didn’t know Euler’s results, you could use the same computation to find Euler’s values: sort the numbers by the sum of fourth powers and use the DIF function to take differences between consecutive elements. Print out the rows for which the difference is zero (that is, that the sum of fourth powers are equal):

/* find the actual pairs of integers that have the same sum */
g = g || f;                    /* concatenate fourth powers onto pairs */
call sort(g, 3);               /* sort by sum of 4th powers */
idx = loc(dif(g[,3])=0);       /* dif=0 <==> consecutive values equal */
print (g[idx-1:idx, ])[c={N1 N2 "Sum"}];

It is fun to think of what someone like Euler might have accomplished if had the tools that we do today!