I tried to repeat this results: 27^5 + 84^5 + 110^5 + 133^5 = 144^5 (Lander & Parkin, 1966)
95800^4 + 217519^4 + 414560^4 = 422481^4 (Roger Frye, 1988, 110 hours this is the only counterexample for values below one million)
The first counterexample can be found by a direct search in several seconds. The second case takes a day to reach 3000 and will likely take years to find the first counterexample. My program iterates over a set of {a1,a2,a3} such that 1<=a1<=a2<=a3. For each triple it calculates the sum of their 4th powers, extracts 4th power root, bring the integer part to a 4th power yielding the nearestPerfectPower. Whilesum is not equal to nearestPerfectPower I choose the next triple {a1,a2,a3} and repeat the steps.

Cases like 1000^4+1^4+1^4<1001^4 slow down the computation and can be removed. I should skip up to 1000^4+1^4+251^4. I might be able to use parity check as well. But even then I don't think I can get the answer fast enough.

How should one efficiently search for a1^4+a2^4+a3^4=a4^4?

  • How is your progress? I made direct search with 300 000 000 comparisons per second, however without any filtering of input numbers.
    – crea7or
    Dec 14, 2016 at 18:10

4 Answers 4


Why not read Frye's article? Here is a digest. You do need modular arithmetic to sieve candidates. Mod 5 works like charm (200 times reduction of the search space). You need to split it in A^4+B^4 and D^4-C^4. The program itself is written in Lisp.

The implementation is reasonably straightforward, while the underlying math is neat. First He sieves out some numbers using modular arithmetics. Then he studies A^4 + B^4 and D^4 - C^4 turning the problem into O(n^2)

First, he notes that mod 5 is very restrictive.

A A^4 A^4 mod 5
0 0 0
1 1 1
2 16 1
3 81 1
4 256 1

Hence, a counterexample must satisfy these two requirements:
A mod 5 = B mod 5 = 0 //A and B are divisible by 5
C^4 mod 5 = D^4 mod 5 = 1 // the reminder of division is 1

Note that D is the largest number and is not divisible by 5, C is not divisible by 5. A < B, A and B are divisible by 5. This reduce the number of possible pairs by a factor of 200. D must be odd using similar considerations for mod 4.

Mod 9 produces only 4 biquadratic yielding 6 congruences of interest. You gain 1.7 improvement from mod 9; 1.4 improvement from mod 13; 1.4 improvement from 29.

If D and C have a common prime factor F then ((AF)^4+(BF)^4 can be eliminated (it was treated earlier). The impact is small, but it helps, nevertheless.

Sieving D^4-C^4 is less efficient, but we still can do certain things. First, we can divide it by 625 before decomposition. Let N be the number to be decomposed. The second constraint is to limit A, B

//I copy-pasted the rest of the article. Don't read it if you want to reinvent the implementation.

(defun  constrained-search  (L  U)    
"Find solution  to  A'  +  B'  +  C '  =  D'    
with  L  <=  D  <  U."       
(decompose-each (constraint-sieve  L  U)))     

This first pass function constraint-sieve iterates on the D variable over the values between L and U which are prime with respect to 10. For each D, it iterates on the C variable over the values less than D which are in proper modulo 625 correspondence with D. When a (D, C) pair passes both the set of modularity tests and the common kctor and prime-factor tests, then the value N = (D' - C')/625 is put into a list. The lists for each choice of D are appended together and returned as the vdue of the function. Here is the code for the first pass function: (defun constraint-sieve (L U) "Apply modular and factor constraints to D and C. Return list of candidates which pass tests."

(loop  for  D  in  (prime-10  L  U)    
(loop for  C  in  (good-e-for-d  D)    
when (and (mod-ok-p  D  C)    
(factor-ok-p  D  C))    
(/  (-  (biquad  D)  (biquad  C))    

The second pass function decompose-each takes the list of candidates generated by the first pass as its argument. It iterates on the variable N over the candidates in the list. For each candidate, it computes the proper limits on A and iterates between them. When it finds an A such that N - A' is a biquadrate, it escapes to the outer level and returns some information about the winner. If a winner is encountered, another function would be needed to fully identify it. Here is the code for the second pass function:

(defun decompose-each (candidates)    
"Attempt  to  decompose each  of  the candidates.   
Return successful decomposition."    
(loop with hd-root  =  (biquad-rt  0.5)    
for  N  in  candidates   
for  a-limit  =  (i-biquad-rt  N)    
for  a-min  =  (ceiling  (*  halflroot a-lim't))    
(loop  for  A from  a-min  below a-limit   
when (biquadrate-p  (-  N  (biquad A)))    
(return-from decompose-each (list  N  A))))))     

How many candidates can we expect the first pass to produce? How many A values does the second pass have to scan for the average candidate? Assume that the candidates are evenly distributed below U', then the average upper limit on A is about 415. (U')'/' = 4U/5 and only 16% of that range is scanned.

  • That's neat. Thank you for the abstract! Our institution doesn't have access to this paper. Where did you get it?
    – sixtytrees
    Jun 29, 2016 at 13:19
  • 1
    @sixtytrees Scientists do research to share their results. They spend time validating hypothesis, writing a manuscript, give it to publisher (and surrender their rights to publisher btw). And publisher says "now pay me money to read it. even if it is your paper. He hides the hard earned knowledge, defeating the whole purpose of the research. Unfair, right? Fortunately, there is a solution. sci-hub.bz (They remove the barriers in the way of science).
    – ygv
    Jun 29, 2016 at 13:29
  • Saved. Great link.
    – sixtytrees
    Jun 29, 2016 at 13:32

I am not sure, how exactly you compute 4th power root, but this operation is very expensive. Likewise computing 4th power is costly. Using a database to precompute 4th powers up to a million might save you a lot of resources.
Likewise, you can start with a perfect power and try subtraction rather than addition. This way you will skip 100^4+1^4+2^4 cases that you mentioned.
Your best bet is to go to a University library and find that paper. (Sorry, I don't have this paper). They must have used some neat tricks that are worth learning about.

  • 4th root is easy: sqrt(sqrt(x)). Likewise, fourth power is easy: x*x*x*x or y*y where y=x*x.
    – user22815
    Jun 29, 2016 at 15:57

Create one data structure that creates all sums of two fourth powers in ascending order, and another one that creates all positive differences of two fourth powers in ascending order. Then compare the numbers from both lists in ascending order. That way you have only O (n^2) numbers to check instead of O (n^3).

Basically, search for a^4 + b^4 = c^4 - d^4. Or let c = d + e, and search for a^4 + b^4 = (d+e)^4 - d^4 = 4d^3e + 6d^2e^2 + 4de^3 + e^3.

To generate a huge list of numbers in sorted order, lookup "priority queue".

Adding a bit about the implementation: Imagine you want to handle sums on the left hand and differences on the right hand up to 10^24. On the left hand side, you have 10^6 ≥ a ≥ b ≥ 0. Use an array for 0 ≤ a ≤ 10^6 storing which b you are using for each a; initially each "a" has b = 0, and you enter the 1,000,001 values a^4 + b^4 into a priority queue (together with each value a). Then you repeatedly take the smallest value from the priority queue, increase the b value in the corresponding array entry by 1, and add the new value a^4 + b^4.

On the right hand side, if you looked at values d^4 - c^4, d could be as large as 63 million. Instead look at values (c + e)^4 - c^4, and you only need e from 1 to 1,000,000.

PS. x^4 = 0 (modulo 16) if x is even, and x^4 = 1 (modulo 16) if x is odd. Therefore in the equation a^4 + b^4 + c^4 = d^4, either a, b, c, d are all even, or d and one of a, b, c are odd. We can ignore the all even solutions, because we get all of those by finding another solution and multiplying all numbers by 2. We can assume that the odd one out of a, b and c is c. So now we look at a^4 + b^4, where a and b are both even, and d^4 - c^4 where d and c are both odd, or (c + e)^4 - c^4, where e is even and c is odd. That saves about 3/4th of all numbers to test.

  • My only concern is that 1M*1M of combinations would require ~1T of entries. Sorting these lists might take a while as well. But it is reasonably easy to implement and O(n^2) is more affordable.
    – sixtytrees
    Jun 27, 2016 at 18:54

Modular arithmetic is the key.

First of all, take the total and add all its bytes together. If any carries occur, add those into the total as well (for instance, F0+F6=E7). Continue until you have a single-byte value. The process works just as well if you start by adding 32-bit words to 32-bit words until you have a 32-bit result and then the 16-bit "top half" to the 16-bit "bottom half" and then the 8-bit "top half" of that result to the 8-bit "bottom half". Since high-level languages don't handle carries very well, use assembler).

This gives you the value of the total, modulo 255 = 3 x 5 x 17. Modulo 3, a 4th power can only be 0 or 1. Modulo 5, a 4th power can only be 0 or 1. Modulo 17, a 4th power can only be 0, 1, 4, 13 or 16. Using these facts you can easily precompute a table which says, for each value between 0 and 255, whether a 4th power (modulo 255) can have that value or not. The table can be a table of 256 bits or 256 bytes: choose whichever is more efficient.

So your algorithm is, given a total, let's say it is 128 bits long:

  1. Add the top 64 bits to the bottom 64 bits, adding in the carry, which will give you a 64-bit number between 0 and 2^64 - 1.
  2. Add the top 32 bits to the bottom 32 bits, adding in the carry, which will give you a 32-bit number between 0 and 2^32 - 1.
  3. Add the top 16 bits to the bottom 16 bits, adding in the carry, which will give you a 16-bit number between 0 and 2^16 - 1.
  4. Add the top 8 bits to the bottom 8 bits, adding in the carry, which will give you an 8-bit number between 0 and 2^8 - 1.
  5. Look up that number in your "might be a 4th power" table, which for convenience has 256 entries for 0 to 255 inclusive (the entry for 255 being equal to the entry for 0).
  6. If the "might be a 4th power" table says that the result might be a 4th power, then test for "4th-powerness" using your existing technique.

Only 4/51, or less than 8%, of your totals will make it as far as step 6, so you have sped up the whole process by a factor of about 12.

As a bonus, with some more precomputation, you can note that the total modulo 255 depends only on the values of your original numbers modulo 255. You can therefore precompute a table of all the 256 x 256 x 256 sets of modulo-255 values of a1, a2, a3 - this is only 16MB or 16Mb, depending how you decide to store it - and use this table to tell you whether it is even worth computing the 4th powers.

You could also try doing things modulo 65535 instead of modulo 255, which removes step 4, leaves only 1/51 = 2% of your totals making it as far as the "fourth root" stage, at the expense of a slightly larger lookup table (64KB or 64Kb).

An even faster filter

Look at the least significant byte only - effectively, the value modulo 256. There are only 18 possible LSBs corresponding to perfect 4th powers, so you reduce your 4th-rooting work to 7% instantly. And this filter can be combined with the modulo-255 or modulo-65535 one. The result will be 181 or 725 times faster than your original enumeration method.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.