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 635,318,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 635318657 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–158 and check to see if the sum of fourth powers contain duplicate values. The number 158 is the greatest integer that is less than the fourth root of 635,318,657.

/* settle conjecture of Euler/Hardy: Is 635,318,657 the smallest
   integer that can be written as the sum of 4th powers in two ways? */
k = floor(635318657##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,561 ordered pairs of integers to consider, and there are 12,560 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!

tags: Just for Fun, Numerical Analysis

3 Comments

  1. Posted October 2, 2012 at 1:27 pm | Permalink

    Fermat's Last Theorem - Go! :D

    • Posted October 2, 2012 at 1:32 pm | Permalink

      Fortunately for the mathematicians, computers still have a hard time with theorems that require checking infinitely many possibilities!

  2. Ian Wakeling
    Posted October 4, 2012 at 4:56 am | Permalink

    For those without IML, something similar to Rick's program can be achieved using PROC SQL. Try the code below to generate the first 16 numbers of this type. Or replace 1000 with 9600 to extend the search to as far as an integer can be exactly represented under Windows. While n=1000 takes less than 1 second, n=9600 takes about 1 min on my machine.

    data n;
      do i=1 to 1000; output;  end;
    run;
    ;
    proc sql;
    create table t as
    select a.i as i1, b.i as i2, i1**4 + i2**4 as s
    from n a, n b where i1 < i2 order by s;
    quit;
    ;
    data _null_;
      set t;
      i3=lag(i1);
      i4=lag(i2);
      if s=lag(s) then put (i3 i4 i1 i2) (5.0) s comma22.0;
    run;

    **************** OUTPUT ******************
       59  158  133  134           635,318,657
        7  239  157  227         3,262,811,042
      256  257  193  292         8,657,437,697
      266  268  118  316        10,165,098,512
      399  402  177  474        51,460,811,217
       14  478  314  454        52,204,976,672
      271  502  298  497        68,899,596,497
      359  514  103  542        86,409,838,577
      512  514  386  584       138,519,003,152
      222  631  503  558       160,961,094,577
      532  536  236  632       162,641,576,192
       21  717  471  681       264,287,694,402
      665  670  295  790       397,074,160,625
      768  771  579  876       701,252,453,457
      354  948  798  804       823,372,979,472
       28  956  628  908       835,279,626,752
    

Post a Comment

Your email is never published nor shared. Required fields are marked *

*
*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <p> <pre lang="" line="" escaped=""> <q cite=""> <strike> <strong>