G. H. Hardy tells the following story about visiting his colleague Ramanujan.

I remember once going to see him when he was ill at Putney. I had ridden in taxi cab number 1729 and remarked that the number seemed to me rather a dull one, and that I hoped it was not an unfavorable omen. “No,” he replied, “it is a very interesting number; it is the smallest number expressible as the sum of two cubes in two different ways.”

This story has become famous, but the rest of the conversation isn’t as well known. Hardy followed up by asking Ramanujan what the corresponding number would be for 4th powers. Ramanujan replied that he did not know, but that such a number must be very large.

Hardy tells this story in his 1937 paper “The Indian Mathematician Ramanujan.” He gives a footnote saying that Euler discovered 635318657 = 158^4 + 59^4 = 134^4 + 133^4 and that this was the smallest number 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?”

Surely someone has since written such a program and settled the question. But as an exercise, imagine the question is still open. Write a program to either come up with a smaller number that answer’s Hardy’s question, or prove that Euler’s number is the smallest one. To make the task more interesting, you might see whether you could do a little pencil-and-paper math up front to reduce the amount searching needed. Also, you might try writing the program in different styles: imperative, functional, etc.

Eureka! I found a smaller number! 635318657

import collections

import math

upper_bound = math.ceil(math.sqrt(math.sqrt(653518657)))

cnt = collections.Counter(x**4+y**4 for x in range(1, upper_bound) for y in range(1, upper_bound) if x >= y)

[x for x in cnt if cnt[x] == 2]

If upper_bound is increased to 1000, the first few results are:

635318657, 3262811042, 8657437697, 10165098512, 51460811217, 52204976672, 68899596497, 86409838577, 138519003152, 160961094577, 162641576192, 264287694402, 397074160625, 701252453457, 823372979472, 835279626752

Here’s my program:

`C=:13 :'#~.0-.~,I."1 x=y'"0 2`

`e4=:3 :0`

n=: }.i.160x

P=: 4 ^~ n

p=:~.,P

A=: +/~p

N=: ~.,A

N#~ 2 < N C A

)

Here’s what it tells me:

`e4''`

635318657

In other words: considering all pairs of fourth powers of integers in the range 1..159, the only sum of pairs that uses more than two of the original integers is 635318657.

That said, note that the original problem assumes we are only considering positive integers. Also, note that 1729 is particularly simple, because 729 is 9^3 and 1728 is 12^3. Note also that once upon a time it memorizing the “times table” for value 1..12 was something everyone did after a few years of school, and memorizing the cubes of 1..12 would be a minor extension of that practice.

P.S. your <code> tag in comments loses indentation.

My script tells me:

635318657 = 59^4 + 158^4, 133^4 + 134^4

3262811042 = 7^4 + 239^4, 157^4 + 227^4

8657437697 = 193^4 + 292^4, 256^4 + 257^4

and so on….

While I’m at it, the smallest number that can be expressed as the sum of three 3rd powers is:

87539319 = 167^3 + 436^3, 228^3 + 423^3, 255^3 + 414^3

You have to love brute force

Very fun. I wonder what Euler would have accomplished if he had access to today’s computational tools? A program and a discussion of the techniques I used are at my blog: Open question in 1937…short SAS program today

Here’s the biggest I found yet:

28^4 + 956^4 = 628^4 + 908^4 = 835 279 626 752

Anyone get their computer to find one over a trillion?

Oh, biggest isn’t the point?

Funny that they both end in 657. Is there a reason for that?

Thanks for the puzzle. I’m learning clojure, so it was a good opportunity!

(defn frth

[x]

(reduce * (repeat 4 x)))

`(defn find-euler-sums`

[]

(let [sums (apply concat (map

(fn [a] (map

(fn [b] (+ (frth a) (frth b)))

(range a 159)))

(range 1 159)))]

(distinct (filter (fn [a] (> (count (filter (fn [b] (= b a)) sums)) 1)) sums))))

searching the space from 1-9000 gives this as the max:

7985644522300177

You can generate arbitrarily big results by scaling a result by n^4. I.e. if

a^4 + b^4 = c^4 + d^4 = S

then, for any integer n,

(na)^4 + (nb)^4 = (nc)^4 + (nd)^4 = S·(n^4)

You could ask, of course, which is the biggest one which is not a scaling of a smaller one.

Has anyone found an example of a number that can be written as a sum of fourth powers in more than two ways? Up to 10^14 there are 69 numbers there that are sums of fourth powers in two ways, but none in 3 or more.

John, @IJ was trying to be funny, so I thin kwe both missed what he was trying to say. There is a typo in your post. The correct number is 635318657, not 653518657.

The last I heard, there were no known examples of numbers that are the sum of two fourth powers in three ways, and it was proven that any example would need to be at least 10^19.

Also, no one knows whether there’s a number that’s the sum of two fifth powers in two ways.

Rick: Thanks. I fixed the error.

@omar I couldn’t calculate it due to floating point arithmetic issues. Perhaps someone using Haskell (fast & arbitrary precision) could set aside an afternoon of computation time to find out

I’m curious what sort of pencil-and-paper calculations John was referring to.

I had my program pre-caculate the fourth powers up to the maximum bound, but I didn’t do anything by hand.

Now that you fixed the error, a new reader (e.g. me) goes back and forth between the post and IJ’s joke, trying to figure out what’s so funny. Making it funnier, in a way.

On OEIS, it’s sequence A016078. The analogous number for fifth powers is not known, but there’s a known lower bound of 2^96.

via @sandeep1chopra: 1729: allowing -ve nos, smallest number’s 91 (6, -5, cubed; 4,3, cubed). Any wonder 91 is divisor of 1729?

For R, I like using the ‘outer’ function for these types of problems. Here’s my ‘quick and dirty’ solution:

`x=1:200`

z=outer(x,x, function(x,y) x^4+y^4)

z=z[upper.tri(z, diag=TRUE)]

names(which(table(z)==2))