Disproving Euler proposition by brute force in C

Question

I wrote a little bit of code to brute force disprove the Euler proposition that states:

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

has no solution when \$a\$, \$b\$, \$c\$, \$d\$ are positive integers.

I'm no mathematician, but reading around this, at least one solution was found by Noam Elkies in 1987 (a = 95800; b = 217519; c = 414560; d = 422481). I wanted to get an idea of how much firepower it would take to solve by brute force, so I wrote the following in C:

#include <stdio.h>
#include <time.h>
#include <math.h>

int prop(long int A, long int B, long int C, long int D) {
    return (pow(A, 4) + pow(B, 4) + pow(C, 4) == pow(D, 4));
}

int main() {
    long int a, b, c, d;
    clock_t t;
    t = clock();

    for (a = 1; a < 100000; a++) {
        for (b = 1; b < 300000; b++) {
            for (c = 1; c < 500000; c++) {
                for (d = 1; d < 500000; d++) {
                    if (prop(a, b, c, d))
                        printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
                }
                if (!(c%1000))
                    printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, ((double)(clock() - t))/CLOCKS_PER_SEC);
            }
            printf("a = %ld, b = %ld, time = %fs\n", a, b, ((double)(clock() - t))/CLOCKS_PER_SEC);
        }
        printf("a = %ld, time = %fs\n", a, ((double)(clock() - t))/CLOCKS_PER_SEC);
    }

}

I ran it for a while, and worked out that it would take roughly \$85 \times 10^6\$ years for it to get to the answer above on my current machine.

I mean, maybe if I waited a few thousand years, I'd have a slightly better machine, but my question is (and I am new to C and computer science in general - so please be gentle); what strategies could I take to make the above code run faster? I thought about threading, and using some sort of bit shifting in place of the pow() calls.

Would there be any way (on my current machine) to get this to run in my lifetime?

Answer 1

Note: At some point, this review drifted into the realm of assembler and GMP. An actual review is at the end of this post, whereas the first section discusses the runtime-problems concerning pow, wrong data types and arbitrary large integers.

No life time for run time

Would there be any way (on my current machine) to get this to run in my lifetime?

There's a great saying in carpentry: measure twice, cut once. It concerns cutting wood or other material, where you have to throw away your resources if you accidentally cut at the wrong place.

A similar saying is there for software engineers: you can't optimize what you can't measure. There are several ways to measure your code, e.g. benchmarking, profiling, or looking at the generated as­sembler to see how many instructions a certain part of your code will take.

Here, we will take the latter route, start with the assembler, take considerations step by step and see where we end up.

A study in assembly

Lets have a look at your code. Well, not yours, but the assembler the compiler generates. You can use gcc -S -O3. On my platform, this results in the following "hot" section in main:

.L6:
        add     rbx, 1
        cmp     rbx, 500000
        je      .L18
.L8:
        mov     rax, QWORD PTR .LC0[rip]
        movsd   xmm0, QWORD PTR [rsp+40]
        movq    xmm1, rax
        call    pow                        ; (1)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+8], xmm0
        movsd   xmm0, QWORD PTR [rsp+48]
        movq    xmm1, rax
        call    pow                        ; (2)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+16], xmm0
        movsd   xmm0, QWORD PTR [rsp+32]
        movq    xmm1, rax
        call    pow                        ; (3)
        mov     rax, QWORD PTR .LC0[rip]
        movsd   QWORD PTR [rsp+24], xmm0
        pxor    xmm0, xmm0
        cvtsi2sd        xmm0, rbx
        movq    xmm1, rax
        call    pow                        ; (4)
        movsd   xmm2, QWORD PTR [rsp+8]
        addsd   xmm2, QWORD PTR [rsp+16]
        movapd  xmm1, xmm0
        movsd   xmm0, QWORD PTR [rsp+24]
        addsd   xmm0, xmm2
        ucomisd xmm0, xmm1
        jp      .L6
        jne     .L6

Even though you might not know assembler, you can see those four calls to pow. The first thing you need to know is that call is slow compared to those other operations. Those four calls happen in the innermost loop. The compiler removed the call to prop and instead replaced it by its code (that's faster).

mov* assigns values to registers, add* adds values, and so on. The registers with xmm* are double precision registers, meant for double variables. So we're basically calling pow with the right values and then add, subtract and modify our small little double values.

Double trouble

But wait a second. We're trying to solve a completely integral problem! Why does our generated program use those registers at all?

This should raise a red flag. And indeed, if we remember pow's signature, it should be clear that it's not the right tool. It takes a double base and exponent, which indicates that it's suitable for terms like \$15.12151^{3.1415926}\$. This is a total overkill for your problem.

Using proper functions

So let's use another pow version instead:

long int pow4(long int x){
    return x * x * x * x;
}

Note that your compiler should create something like this from that:

movq    %rdi, %rax
imulq   %rdi, %rax
imulq   %rax, %rax
ret

but if your compiler doesn't recognize this potential (micro) optimization, you can use

long int pow4(long int x){
    const long int t = x * x;
    return t * t;
}

instead.

We also need to change prop:

int prop(long int A, long int B, long int C, long int D) {
  return (pow4(A) + pow4(B) + pow4(C) == pow4(D));
}

Allright. Now, before I show the times of the new program, let's check the output of your old one:

a = 1, b = 1, c = 1000, time = 114.156000s

That's when I hit ^C. How does the one using pow4 hold up?

a = 1, b = 1, c = 1000, time = 0.296000s
a = 1, b = 1, c = 2000, time = 0.578000s
a = 1, b = 1, c = 3000, time = 0.859000s
a = 1, b = 1, c = 4000, time = 1.140000s
a = 1, b = 1, c = 5000, time = 1.421000s
a = 1, b = 1, c = 6000, time = 1.703000s
a = 1, b = 1, c = 7000, time = 1.984000s
a = 1, b = 1, c = 8000, time = 2.265000s
a = 1, b = 1, c = 9000, time = 2.546000s
a = 1, b = 1, c = 10000, time = 2.828000s
a = 1, b = 1, c = 11000, time = 3.109000s
a = 1, b = 1, c = 12000, time = 3.390000s
a = 1, b = 1, c = 13000, time = 3.687000s
a = 1, b = 1, c = 14000, time = 3.968000s
a = 1, b = 1, c = 15000, time = 4.250000s
a = 1, b = 1, c = 16000, time = 4.531000s

Which is 0,2% of your original time, or a 500x speedup.

However, this comes at a cost: pow4(500000) is too large for a int64_t, since \$\log_2(500000^4) \approx 76\$. The greatest number you could check with a uint64_t is 65535, \$2^{16}-1\$, which shouldn't be very surprising. As the standard does not provide int128_t or similar, you should make sure that your numbers don't exceed those bounds.

You can either write your own large integer logic for this, or use GMP.

Proper bounds and parameter estimation

Next up, you can increase the lower bounds of b and c, so that \$a \le b \le c\$. And for d, well, if we have a, b, c, then there is only one solution for d. We can directly search for that solution with binary search.

The binary search makes a \$\mathcal O (n^3 \log n)\$ algorithm from your current \$\mathcal O (n^4)\$ one, which should provide a lot more speed than the previous speedup.

Even better, if you used the appropriate bounds for a, b and c, we can bound d by

$$d^4 = a^4 + b^4 + c^4 \le c^4 + c^4 + c^4 = 3c^4$$

and therefore get

$$c \le d \le \sqrt[4]{3}\,c.$$

With the proper binary algorithm,you can finish the first a=1,b=1 case quickly:

…
a = 1, b = 1, c = 481000, time = 0.031000s
a = 1, b = 1, c = 482000, time = 0.031000s
a = 1, b = 1, c = 483000, time = 0.031000s
a = 1, b = 1, c = 484000, time = 0.031000s
a = 1, b = 1, c = 485000, time = 0.031000s
a = 1, b = 1, c = 486000, time = 0.031000s
a = 1, b = 1, c = 487000, time = 0.031000s
a = 1, b = 1, c = 488000, time = 0.031000s
a = 1, b = 1, c = 489000, time = 0.031000s
a = 1, b = 1, c = 490000, time = 0.031000s
a = 1, b = 1, c = 491000, time = 0.031000s
a = 1, b = 1, c = 492000, time = 0.031000s
a = 1, b = 1, c = 493000, time = 0.031000s
a = 1, b = 1, c = 494000, time = 0.031000s
a = 1, b = 1, c = 495000, time = 0.031000s
a = 1, b = 1, c = 496000, time = 0.031000s
a = 1, b = 1, c = 497000, time = 0.031000s
a = 1, b = 1, c = 498000, time = 0.031000s
a = 1, b = 1, c = 499000, time = 0.031000s
a = 1, b = 1, time = 0.031000s

Which brings us back into the realm of your lifetime.

Exercise

Write a function, that given a, b and c checks whether there exist a d, such that your property holds. It should return -1 if there does not exist such a d, and the d otherwise.

Use that function in your code. Make sure that you need roughly \$\log d_{\text{max}}\$ iterations in that function.

Important remark about integer sizes

Keep in mind that long int is usually just a 64 bit integer, which means that the largest integer you can store is \$2^{63}-1\$. Integer types with more bits have greater bounds, but are platform specific. Also, multiplication can be a tad slower, since multiplying 128bit numbers isn't as easy as multiplying 64bit numbers.

See the next section how to get multiplications down.

An actual review

Our pow4 is now essentially two multiplications. However, we're still using pow4 too often. After all, we don't need to recalculate \$a^4\$ in every iteration. The compiler happily does, since it doesn't optimize aggressively enough.

Which brings us to the actual review: your code is cleanly written, easy to read and to understand. Unfortunately, well-written, modular code often doesn't squeeze the last bit (heh) out of your hardware, unless your compiler/runtime is very smart (and thus often expensive).

So let's get back to the drawing board for a final review of your code:

Includes

#include <stdio.h>
#include <time.h>
#include <math.h>

I would sort them by name, but that's fine. You don't include anything that's not necessary, nor did you forget something (and got saved by a non-standard compliant compiler).

Declaration

int main() {
    long int a, b, c, d;
    clock_t t;
    t = clock();

Depending on whether you write ANSI-C or C99, I would defer declaration of variables as long as possible. For example, at the moment it's easy to accidentally change c to some bogus value, or forget a { and accidentally check the prop after the loops or similar:

for (a = 1; a < 100000; a++) 
    for (b = 1; b < 300000; b++) 
        for (c = 1; c < 500000; c++) 
            for (d = 1; d < 500000; d++) 
                printf("inner loop");
                if (prop(a,b,c,d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);

Whoops. The if doesn't get checked, and you don't get a warning (in older compiler versions; new ones do warn about possible whitespace issues). If you declare your variables later (e.g. C99-style), errors like that cannot happen (although it introduces possible shadowing):

for (long int a = 1; a < 100000; a++) 
    for (long int b = 1; b < 300000; b++) 
        for (long int c = 1; c < 500000; c++) 
            for (long int d = 1; d < 500000; d++) 
                printf("inner loop");
                if (prop(a,b,c,d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd %ld\n", a, b, c, d);

This will now lead to a compiler error, since a, b and so on are out of scope. Either way, that depends on the language standard you want to use. Some people prefer one way, others the other one. Choose yours.

Types

Given that all values should be strictly greater than zero, long int is not the appropriate type, as it can be negative. We should accommodate that. However, instead of using long unsigned int through­out our code, let's use a type synonym in case we want to change it later to a type with a greater range:

 typedef long unsigned int Number;

You can probably come up with a better name.

Cache results (by hand)

One thing that strikes me most is that you recalculate \$a^4\$ and so on every time. We can easily treat this with more variables (using your declaration style):

int main() {
    long int a, b, c, d;
    long int a4, b4, c4, d4; // new variables
    clock_t t;
    t = clock();

    for (a = 1; a < 100000; a++) {
        a4 = pow4(a);                             // remember
        for (b = a; b < 300000; b++) {
            b4 = pow4(b);                         // remember
            for (c = b; c < 500000; c++) {
                c4 = pow4(c);                     // the fourth power
                for (d = c; d < 500000; d++) {
                    d4 = pow4(d);                 // of this member
                    if (a4 + b4 + c4 == d4)
                        printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
…

Remember how I said that nicely written, modular code isn't often optimal? This is one of those un­fortunate examples where you have to help the compiler (unless you know exactly what optimization flags you have to use or your compiler is overly aggressive). The prop is gone, the calls to pow4 are now in your loop.

But the compiler cannot make a mistake here anymore: it's now very clear that a4 doesn't need to be recalculated 300000*500000*500000 times.

That being said, we should apply the other suggestions like the type synonym and the late declaration:

typedef long unsigned int Number;

int main() {
    clock_t t;
    t = clock();

    for (Number a = 1; a < 100000; a++) {
        const Number a4 = pow4(a);                         // remember
        for (Number b = a; b < 300000; b++) {
            const Number b4 = pow4(b);                     // remember
            for (Number c = b; c < 500000; c++) {
                const Number c4 = pow4(c);                 // the fourth power
                for (Number d = c; d < 500000; d++) {
                    const Number d4 = pow4(d);             // of this member
                    if (a4 + b4 + c4 == d4)
                        printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
…

While const isn't necessary here, it will make sure that we don't change our cached results accidentally.

The time has come

Although our code is now more verbose, there is one small piece of code that repeats itself three times throughout your main:

((double)(clock() - t))/CLOCKS_PER_SEC)

That's quite hard to read, isn't it? It's a perfect candidate for a function:

static inline seconds_since(clock_t t){
    return ((double)(clock() - t))/CLOCKS_PER_SEC;
}

This changes your printf from

printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, ((double)(clock() - t))/CLOCKS_PER_SEC);

to

printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, seconds_since(t));

Ah. Much easier to read. That's what inline functions are for. Note that any sophisticated compiler should inline that function anyway, so you may also drop inline if you don't want to use C99.

Answer 2

Well let's start with this:

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {

let's ignore d for now. What do you do here? You check 1, 1, 1, then 1, 1, 2, then 1, 1, 3, ... up to 1, 1, 499999. Then you start over at 1, 2, 1. But you already checked 1, 1, 2, so why are you checking 1, 2, 1? You could go straight to 1, 2, 2. That doesn't save you much for these low numbers, but believe me, when you get to big numbers it adds up.

In short: a, b, c should be nondecreasing. We can achieve that by starting b at a, and starting c at b, so b is never smaller than a, and c is never smaller than b. So immediately you can get rid of about half the work you're doing with

for (a = 1; a < 100000; a++) {
    for (b = a; b < 300000; b++) {
        for (c = b; c < 500000; c++) {

Next, consider d.

Of course we can immediately make the same optimization we just made for a, b and c, since d will never be smaller than any of them, and c is always the largest.

Also, once d⁴ is larger than the sum, we can stop incrementing d because it's only going to get bigger.

So that will save a lot of time right there. But we can do way better than that.

The question you are asking is "do these four numbers have the sum property?" but the question you should be asking is "does a⁴ + b⁴ + c⁴ equal any fourth power?" If it does, then you can easily compute d much faster than trying all possible fourth powers. So, can you write a fast method that tells you if a particular sum is a fourth power or not?

If you know how to take a square root, what you can do is take the sum, take the square root twice, and then square the result twice. If you get back the original sum, then it was a fourth power, and if you don't, then it wasn't.

If you don't know how to take a square root, you can do the following logic: we have sum; is 128⁴ equal to sum? No. Is it bigger? Yes. Then next try 64⁴. Is it equal to sum? No. Is it smaller? Yes. Then try 96⁴, and so on. Binary search for the result until you find it, or until there are no more numbers in your range to check. That does only a tiny handful of computations, compared to trying thousands upon thousands of possible fourth powers.

Answer 3

All of the answers given (when this 'answer' was first written) ignore an important constraint: native types. In C, a long int is a 32-bit (or greater) signed type, meaning the largest positive value that can be (counted upon to be) represented is \$2^{31}-1\$. The largest possible input to a function which calculates a fourth-power and [[can be counted on to]] return a correct long int would be \$\left\lfloor\sqrt[4]{2^{31}-1}\right\rfloor=215\$. Even if you move to unsigned long long int, a 64-bit (or greater) unsigned type, \$\left\lfloor\sqrt[4]{2^{64}-1}\right\rfloor=65\,535\$. Floating-point types wouldn't get you any further: only the mantissa bits are of value because of the precision needed. (If you have access to a 128-bit type, that would suffice, as \$\left\lfloor\sqrt[4]{2^{128}-1}\right\rfloor=4\,294\,967\,295\$.)

If native types won't fully suffice, you will need a 'big number library' that handles such values. A word of warning: having to call out to a library will most likely incur a performance hit.

All that said, it may still be possible to write the code mostly with native types (under 128 bits) in such a way as to add and compare with only the lower bits, but you would have to be very careful, and you have the possibility of false-positives. At that point you could call in the 'big number library' to verify the find.

Answer 4

Some interesting things no one have mentioned about this, but it is an improvement when looking for the smallest solution:

If the \$\gcd(a,b,c) \neq 1\$ then \$a^4\$,\$b^4\$,\$c^4\$ and \$d^4\$ are all divisible by that gcd to the 4th power, giving a smaller solution.

Therefore, at least one of \$a\$,\$b\$, and \$c\$ must be odd, since if they are all even they share a common divisor of 2, meaning the \$\gcd(a,b,c) \geq 2\$. By checking that they aren't all even, before performing the loop for d, you will cut your possible search set by a large amount.

Along these same lines, if \$d\$ is even, then 2 must divide \$a^4 + b^4 + c^4\$, but we already mentioned that they can't all be even to be the smallest solution. If two (or zero) out of the three variables are even then the sum of their 4th powers is going to be odd (even + even + odd = odd; odd + odd + odd = odd), thus if d is even then 1 number must be even and the other 2 are odd.

If \$d\$ is odd, then 2 cannot divide \$a^4 + b^4 + c^4\$. If only 1 of the numbers is even, then then the sum of their 4th powers is going to be even (even + odd + odd = even). By this either they are all odd or only 1 is odd.

Using this you can decide if D is odd or even based on the set of \$a,b,c\$ cutting your solution set in half by tailoring the innermost loop accordingly. In conjunction with the increasing values d>c>b>a, this should improve your processing by a large margin.

Answer 5

I'm not able to comment on vnp's solution, but vnp was right the first time: you can brute force it in \$O(n^2\log(n))\$ time and \$O(n)\$ space. You don't need \$O(n^2)\$ space because you don't have to store the whole list of \$a^4+b^4\$ or \$d^4-c^4\$ upfront. Instead you only need to be able to list the values of \$a^4+b^4\$ and \$d^4-c^4\$ in ascending order and then match up the two values.

To list the values of \$a^4+b^4\$, you can maintain, for each \$a\$, the smallest value of \$b\$, call it \$b_a\$, such that \$a^4+b_a^4\ge k\$, where \$k\$ is the previous value of \$a^4+b^4\$ that you read off. To read off the next value, you need to find the value of \$a\$ such that \$a^4+b_a^4\$ is the smallest, which you can do in \$O(\log(n))\$ time using a heap or something. Then you increment \$b_a\$ and you are ready to read off the next value of \$a^4+b^4\$.

This is an example C++11 program (using a GNU extension to get 128 bit integers) that finds the solution in about 7 hours on my PC. (This is just to illustrate the method - it is not very efficient in other ways.)

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <queue>
#include <vector>
using std::priority_queue;
using std::vector;

//typedef long long int bigint;
typedef __int128 bigint;

// Look for a^4+b^4+c^4 = d^4 by comparing
// ascending list of a^4+b^4 with ascending list of d^4-c^4

vector<bigint> list;
vector<int> diffptr,sumptr;

bool sumcmp(int a0,int a1){
  int b0=sumptr[a0];
  int b1=sumptr[a1];
  return list[a0]+list[b0]>list[a1]+list[b1];
}

bool diffcmp(int c0,int c1){
  int d0=diffptr[c0];
  int d1=diffptr[c1];
  return list[d0]-list[c0]>list[d1]-list[c1];
}

int main(int ac,char**av){
  int a,c,i,n;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint sumval,diffval;

  list.resize(n);
  diffptr.resize(n);
  sumptr.resize(n);

  priority_queue<int,vector<int>,decltype(&sumcmp)> sumnext(sumcmp);
  priority_queue<int,vector<int>,decltype(&diffcmp)> diffnext(diffcmp);

  for(i=0;i<n;i++)list[i]=bigint(i)*i*i*i;
  for(a=1;a<n;a++){sumptr[a]=a;sumnext.push(a);}
  for(c=0;c<n-1;c++){diffptr[c]=c+1;diffnext.push(c);}

  while(!sumnext.empty()&&!diffnext.empty()){
    a=sumnext.top();sumval=list[sumptr[a]]+list[a];
    c=diffnext.top();diffval=list[diffptr[c]]-list[c];
    if(sumval==diffval){printf("BINGO %d^4+%d^4+%d^4=%d^4\n",a,sumptr[a],c,diffptr[c]);fflush(stdout);}
    if(sumval<=diffval){
      sumnext.pop();
      sumptr[a]++;if(sumptr[a]<n)sumnext.push(a);
    }else{
      diffnext.pop();
      diffptr[c]++;if(diffptr[c]<n)diffnext.push(c);
    }
  }

}

Answer 6

Read the source

Since you've already identified the author, why not read the original paper?

Some misconceptions

In fact, the numbers found by Elkies were not the ones cited in the question. They were found by Roger Frye using a computer. How did he do that back in 1988? It's mentioned in the paper:

Postscript.

While our first counterexample $$(A, B, C;D) = (2682440,15365639,18796760; 20615673)$$ to Euler's conjecture still seems beyond the range of reasonable exhaustive computer search, there remained the possibility that smaller solutions may be found by such a search. Shortly after hearing of the first solution, Roger Frye of Thinking Machines Corporation asked whether it was minimal; I did not know, but suggested how one might exhaustively search for smaller solutions: eliminating common factors and permuting \$A\$, \$B\$, \$C\$ if necessary, we may take \$D\$ odd and not divisible by 5, and \$C < D\$ such that \$D^4 - C^4\$ is divisible by 625 and satisfies several other congruence and divisibility properties, and for each such \$D\$ and \$C\$ look for a representation of \$D^4 — C^4\$ as \$A^4 + B^4\$ with \$A\$, \$B\$ divisible by 5. Frye translated this into a computer program and ran it on various Connection Machines for about 100 hours to find the minimal counterexample to Euler's conjecture: $$95800^4 + 217519^4 + 414560^4 = 422481^4.$$ He continued the search and found that this solution is unique in the range \$D < 10^6\$. This solution appears on the parameterization (6) with \$(m, n) = (20,-9)\$. We include Frye's result with his permission.

Source: Elkies, Noam D. "On \$A^4+B^4+C^4=D^4\$." Mathematics of Computation (1988): 825-835.

For those unfamiliar with the late great Thinking Machines Corporation, it was a company in the 1980s and 1990s that produced massively parallel computers.

Taken together, these hints suggest that a multithreaded approach using an exhaustive search with these number theoretic tests should definitely be doable on modern machines in the range of hours rather than years.

Answer 7

I recommend to rewrite a problem as finding \$a,b,c,d\$ such that \$a^4 + b^4 = d^4 - c^4\$. Now you may only operate in pairs of powers. Building a table of sums takes \$O(n^2)\$. Building the table of differences also takes \$O(n^2)\$. Sorting them takes \$O(n^2\log{n})\$. Matching sorted tables takes \$O(n^2)\$ in a linear merge-like scan.

Overall complexity is obviously \$O(n^2\log{n})\$.

Answer 8

This is a different method from the other (priority queue) answer I submitted earlier. It is simpler and achieves \$O(n^2\log(n))\$ time and \$O(n)\$ memory with a better constant factor. You choose some intervals \$[A_i,A_{i+1})\$ covering \$[0,n^4)\$ and then separately (for each \$i\$) find \$a,b\$ such that \$a^4+b^4\in[A_i,A_{i+1})\$ and \$c,d\$ such that \$d^4-c^4\in[A_i,A_{i+1})\$. Then you sort the two lists of \$a^4+b^4\$ and \$d^4-c^4\$, and check them for an intersection.

To make sure that it runs in the right time order, you need to make sure that the number of \$a,b\$ with \$a^4+b^4\in[A_i,A_{i+1})\$ is proportional to \$n\$ and roughly independent of \$i\$. You can do this by using \$A_i=C\cdot n^2\cdot i^2\$ for some constant \$C\$. The optimal value of \$C\$ will depend on your specific hardware characteristics.

If you do this (first program below) then the running time to find the above solution \$95800^4+414560^4+217519^4=422481^4\$ seems to be a little over an hour (on my desktop computer). This uses __int128, a GCC extension giving 128 bit integers.

For fun, I've also modified it to use the simplest (and most useful) modulo relations: \$a,b\$ are multiples of \$5\$, and \$d^4-c^4\$ is a multiple of \$625\$. This is the second program below and it finds the solution in about four minutes.

I've also modified it to be multi-threaded. This is straightforward, because the chunks \$[A_i,A_{i+1})\$ can be treated separately. This is the third program below and it finds the solution (on my hex core machine) in about 32 seconds.

Another technique that speeds things up is to reduce the elements of the two lists modulo \$2^{64}\$ before sorting and comparing them. If there is a match between \$a^4+b^4\$ and \$d^4-c^4\$ mod \$2^{64}\$ then you need an extra check that it is genuine, meaning it arises from an unreduced match. (It is likely to be genuine, so the extra check takes negligible extra time.) Because most of the time is spent sorting, it is better to sort 64 bit quantities than 128 bit quantities. This reduces the memory requirement by a factor of (nearly) 2 and the running time by about 30%, making it about 23 seconds to the first solution on a hex core machine. This is the fourth program below.

(It does make assumptions like casting the 128 bit numbers to 64 bits is the same as taking the bottom 64 bits, and long long int is 64 bits, and I've left some Unix-specialised timing code, so it's not all that pretty or portable.)

Additional notes:

• Programs 2 and 3 edited to incorporate Edward's improvement.

• If you run program 4 with n=3000000, it finds the new solution \$673865^4+1390400^4+2767624^4=2813001^4\$ after about 5 minutes (on a hex core machine). If you run it with n=9000000 and MAXMEM=12e9, it finds \$1705575^4+5507880^4+8332208^4=8707481^4\$ after about 80 minutes.

• It may be possible to enhance this interval method to run in \$O(n^2)\$ time and \$O(n^{2/3})\$ memory by using continued fractions. This is quite fiddly (meaning annoying constant factor overheads) and would probably only be useful for \$n\ge10^8\$ or so, assuming memory of the order of 32GB.

• The state of the art for this problem, devised by Robert Gerbicz, is recorded here and here (Russian). Their program uses many more "modulo filters", and uses the Chinese remainder theorem instead of meet-in-the middle, which means it uses much less memory as well as being faster. They have found solutions up to \$d=1871713857\$.

Program 1 - basic chunking

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4).

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

vector<bigint> pow4;
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,s0,s1;
  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=bigint(min(5.*n,MAXMEM/26.));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){
    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)

    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=1;a<n&&pow4[a]<ub;a++){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=min(int(sqrt(sqrt(ub-pow4[a]))),n);
      for(b=bmin;b<bmax;b++)list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    for(c=0;c<n-1;c++){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(d=dmin;d<dmax;d++)list1.push_back(pow4[d]-pow4[c]);
    }
    std::sort(list1.begin(),list1.end());

    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }
  }
}

Program 2 - as program 1, but also use modulo 5 cutdown

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

#define MAXMEM 4e9 // Allow the use of up to (approx) this memory in bytes

vector<bigint> pow4;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,r,s0,s1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  pow4.resize(n);

  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint((MAXMEM-2*n*sizeof(int))*2.8));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));
  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(chunk/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

Program 3 - as program 2, but also use multithreading

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.
// Parallelise over [A_i,B_i).

// Compile with -pthread flag

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
#include <thread>
#include <mutex>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

int n;// Maximum value of d
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
const int nthreads=12;
bigint chunk,chunksize;
vector<bigint> pow4;
std::mutex chunkserver,printer;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  std::lock_guard<std::mutex> lock(printer);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

bigint nextchunk(){// Returns next chunk, or -1 if none left
  std::lock_guard<std::mutex> lock(chunkserver);
  if(chunk>=bigint(n)*n)return -1;
  bigint chunkcopy=chunk;
  chunk+=chunksize;
  return chunkcopy;
}

void searcher(){
  int a,b,c,d,i,j,r,s0,s1;
  bigint ch;
  vector<bigint> list0,list1;
  while((ch=nextchunk())>=0){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=ch*ch;
    bigint ub=min(ch+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(ch/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

int main(int ac,char**av){
  int i,t;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint(MAXMEM*2.8/nthreads));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  vector<std::thread> thr(nthreads);
  for(t=0;t<nthreads;t++)thr[t]=std::thread(searcher);
  for(t=0;t<nthreads;t++)thr[t].join();
}

Program 4 - as program 3, but reduce 128-bit integers to 64 bits

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.
// Parallelise over [A_i,B_i).

// http://codereview.stackexchange.com/questions/145221/disproving-euler-proposition-by-brute-force-in-c

// Compile with -pthread flag

#include <cstdio>
#include <cmath>
#include <ctime>
#include <sys/time.h>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <vector>
#include <thread>
#include <mutex>
using std::vector;
using std::min;
using std::max;

typedef __int128 bigint;// GCC extension: 128 bit integer
typedef unsigned long long int medint;// Need this to be 64 bits

int n;// Maximum value of d
double tim0,tim1,dtim;
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
const int nthreads=12;
bigint chunk,chunksize;
vector<bigint> pow4;
std::mutex chunkserver,printer;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}
double cpu(){return clock()/double(CLOCKS_PER_SEC);}
double realtime(){
  struct timeval tv;
  gettimeofday(&tv,0);
  return tv.tv_sec+tv.tv_usec/1e6;
}

void printdecode(bigint cv){// Find a,b,c,d given a possible common value, cv=a^4+b^4=d^4-c^4
  int a,b=0,c,d=0;// Initialise to prevent unnecessary compiler warnings
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  if(a==n)return;
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(d>=n)return;
    if(pow4[d]-pow4[c]==cv)break;
  }
  if(c==n-1)return;
  std::lock_guard<std::mutex> lock(printer);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4 @ CPU=%.2fs Real=%.2fs\n",a,b,c,d,cpu(),realtime()-tim0);fflush(stdout);
}

bigint nextchunk(){// Returns next chunk, or -1 if none left
  std::lock_guard<std::mutex> lock(chunkserver);
  if(chunk>=bigint(n)*n)return -1;
  bigint chunkcopy=chunk;
  chunk+=chunksize;
  {
    double tim2=realtime()-tim0;
    if(tim2>=tim1){
      tim1=tim2+dtim;dtim=min(1.2*dtim,100.);
      double adjchunk=chunk-(nthreads+1)/2.*chunksize;// Estimated chunkage completed
      double p=adjchunk/double(bigint(n)*n);// Approx proportion complete (not strictly correct)
      printf("%11.2f   %8.3f%%   %12g   %11.2f  %11.2f\n",tim2,p*100,sqrt(max(adjchunk,0.)),tim2/p,tim2/p*(1-p));fflush(stdout);
    }
  }
  return chunkcopy;
}

void searcher(){
  int a,b,c,d,i,j,r,s0,s1;
  bigint ch;
  vector<medint> list0,list1;
  while((ch=nextchunk())>=0){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=ch*ch;
    bigint ub=min(ch+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=max(int(sqrt(sqrt(lb+pow4[c]))),c+1);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(ch/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){// Match modulo 2^64: find full 128-bit value to see if this is a genuine match
        for(a=5;a<n&&pow4[a]<ub;a+=5){// Ugly code repeat
          int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
          int bmax=int(sqrt(sqrt(ub-pow4[a])));
          r=(-bmin)%5;if(r<0)r+=5;
          for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1)&&medint(pow4[a]+pow4[b])==list0[i])printdecode(pow4[a]+pow4[b]);
        }
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

int main(int ac,char**av){
  int i,t;

  n=(ac>1?atoi(av[1]):500000);
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}
  assert(sizeof(medint)>=8);

  tim0=realtime();
  tim1=0;dtim=1;
  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(600)*n,bigint(MAXMEM*5.8/nthreads));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  vector<std::thread> thr(nthreads);
  for(t=0;t<nthreads;t++)thr[t]=std::thread(searcher);
  for(t=0;t<nthreads;t++)thr[t].join();
  printf("Total CPU %.2fs\n",cpu());
  printf("Total real %.2fs\n",realtime()-tim0);
}

Answer 9

Multiplication is not needed.

Consider how to calculate a series of squares 1,4,9,16

unsigned sq = 1;
for (unsigned x=1;; x++) {
  sq += x + x + 1;
  printf("%u\n", sq);
}

Consider how to calculate a series of cubes 1,8,27,64

unsigned sq = 1;
unsigned cube = 1;
for (unsigned x=1;; x++) {
  cube += (sq + sq + sq) + (x + x + x) + 1;
  sq += x + x + 1;
  printf("%u\n", cube);
}

A similar operation for 4th power.

---

Wider, but not arbitrarily wider math needed

To reach a = 95800; b = 217519; c = 414560; d = 422481, 80 bit math is needed. Suggest using a complier that supports uint128_t or equivalent.

---

Just integer addition and compare are needed

I would incrementally calculate a**4 + b**4 + c**4 and d**4, increment the smaller of the two. Of course, when they are equal - voila!

---

Sample code. Use factor = 1 for the final code. Use factor >= 7 for testing.

Using a factor allows code to run just a portion of the values needed for testing. A factor of 50 took a few minutes. A factor of 20 took a 4 hours. Expect a factor of 1 would take years (5?). When factor < 7, you need uint128t type that is really 128 bits, unlike the uintmax_t place-holder used here.

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

// Adjust to some 128 bit type
// The below is a place-holder for a true 128-bit type
typedef uintmax_t uint_fast128t;

typedef struct {
  uint_fast128t x4;  // 4th power of the value
  uint_fast128t x3;  // cube of the value
  uint_fast64_t x2;  // square of the value
  uint_fast32_t x1;  // location of OP's value a,b,c or d
                     // Assume a,b,c,d are 20 bits or less.
}quad;

#define quad_next(x) { \
  x.x4 += (uint_fast128t) x.x3*4 + x.x2*6 + 4*x.x1 + 1; \
  x.x3 += (uint_fast128t) x.x2*3 + x.x1*3 + 1; \
  x.x2 += (uint_fast64_t) x.x1*2 + 1; \
  x.x1++; \
  }

#define factor 20
uint_fast32_t a_max = 95800 / factor;
uint_fast32_t b_max = 217519 / factor;
uint_fast32_t c_max = 414560 / factor;
uint_fast32_t d_max = 422481 / factor;

void Elkies_B(quad a, quad b) {
  quad c = b;
  quad d = c;
  c.x4 += b.x4;

  // Could use a more efficient step here
  while (d.x4 < c.x4) {
    quad_next(d);
  }
  while (c.x1 <= c_max) {
    while (d.x4 < c.x4) {
      quad_next(d);
    }
    if (d.x4 == c.x4) {
      // Success! Print a,b,c,d
      printf("+ a:%lu b:%lu c:%lu d:%lu\n", //
              (unsigned long) a.x1, (unsigned long) b.x1, //
              (unsigned long) c.x1, (unsigned long) d.x1);
      exit(0);
    }
    quad_next(c);
  }
}

void Elkies_A(quad a) {
  quad b = a;
  b.x4 += a.x4;
  while (b.x1 <= b_max) {
    Elkies_B(a, b);
    quad_next(b);
  }
}

void Elkies(void) {
  quad a = { 1, 1, 1, 1 };
  while (a.x1 <= a_max) {
    printf("- a:%lu %llu %ju %ju\n", //
            (unsigned long) a.x1, (unsigned long long) a.x2, //
            a.x3, a.x4);
    fflush( stdout);
    Elkies_A(a);
    quad_next(a);
  }
  printf("Failed\n");
}

int main() {
  Elkies();
}

Answer 10

The solution you are talking about was not found by Elkies but by Frye. This is described in [Elkies]:

Euler conjectured in 1769 that the Diophantine equation \$A^4 + B^4 + C^4 = D^4\$,(...)has no solution in positive integers.

(...) eliminating common factors and and permuting \$A,B,C\$ if necessary, we may take \$D\$ odd and not divisible by \$5\$ and \$C \lt D\$ such that \$D^4-C^4\$ is divisible by \$625\$ and satisfies several other congruence and divisibility properties, and for each such D and C look for a representation of D^4-C^4 as \$A^4+B^4\$ with \$A,B\$ divisible by \$5\$. Frye translated this into a computer program and ran it (...) for about 100 hours [in 1988] to find the minimal counterexample to Euler's conjecture

Roger Frye published a detailed description of his algorithm here:[Frye].

The first investigation (for values less than 10000) was made by [Ward]. He didn't use a computer. The paper isn't available for free. But the first page containing the important theorems one and two can be found in Duke Mathematical Journal and a summary that is referenced from Гипотеза Эйлера can be found here.

Another investigation (now by computer) for values less than 220 000 was published here: [Lander,Parkin,Selfridge]

More recent investigations: [Ekl] and [Bernstein].

The most advanced search (values smaller 2 000 000) I found is published by Robert Gerbicz' page: it is the C-program euler413.c.

_{Note from @AlexSelby: To make the code work on a 64-bit system you need to change the five lines `multipliers17=(unsigned int**) (malloc) (17*17*sizeof(unsigned int));` and so on, replacing `sizeof(unsigned int)` with `sizeof(unsigned int*)` (and similarly for the following four lines).}

---

[Bernstein], Bernstein, D.J.; ENUMERATING SOLUTIONS TO p(a) + q(b) = r(c) + s(d); MATHEMATICS OF COMPUTATION, Vol 70, N0 233, pp 389-394

[Ekl], Ekl, R.L.; NEW RESULTS IN EQUAL SUMS OF LIKE POWERS; MATHEMATICS OF COMPUTATION, Vol 67, No 223, July 1998, pp 1309-1315

[Elkies] Elkies, Noam D.; On \$A^4 + B^4 + C^4 = D^4\$; MATHEMATICS OF COMPUTATION, 51/184, OCT. 1988, 825-835

[Frye] Frye, Roger; Finding $95800^4 + 217519^4 + 414560^4 = 422481^4$ on the Connection Machine

[Lander,Parkin,Selfridge] Lander,L. J.; Parkin T. R.; Selfridge, J. L.; A Survey of Equal Sums of Like Powers; Math. Of Computation, Vol 21, 1967

[Ward] Ward, Morgan. Euler’s problem on sums of three fourth powers. Duke Math. J. 15 (1948), no. 3, 827--837
ward.txt is a summary and here is the first page of Ward's article

---

---

Just for fun!

A compiled version for Windows can also be found on Robert Gerbicz' page: euler413.exe. I tried it out. I know nothing about the program, its parameters and its output. I got the following:

C:> euler413.exe
I haven't found unfinished work!
Please give R parameter, the Range will be R*10240000: 10
Do you want to start an exhaustive search in this Range or
only to search for special solutions?
( y=full search, n=special search ) n
Give the start value of a0. It should be 0<=start_a0<=16384: 0
Give the end value of a0. It should be start_a0<=end_a0<=16384: 16384
Building up some tables
Done
Testing: a0=1
Finished: a0=1,Range=102400000,type=0,Time: 0h0m0s,Date: Fri Nov 04 18:20:33 2016
Testing: a0=9
Finished: a0=9,Range=102400000,type=0,Time: 0h0m1s,Date: Fri Nov 04 18:20:34 2016
....

The program creates a lot of output to the screen and the following files:

euler413work.txt
stat_euler(4,3,1).txt
results_euler(4,1,3).txt

euler413work.txt is a status file. If you stop the program with <Ctrl+C> and start it later in the same directory it will continue the work at the point you stopped it. stat_euler(4,3,1).txt records the timing output that you can also see on the screen. In results_euler(4,1,3).txt the solution founded so far are stored.

After about 9 minutes the results-file contained the following lines:

Solution found! 98438073^4=50681927^4+96592480^4+22321400^4
Solution found! 17321721^4=8918279^4+16996960^4+3927800^4
Solution found! 40980657^4=21099343^4+40212320^4+9292600^4
Solution found! 20615673^4=15365639^4+18796760^4+2682440^4  (primitive)
Solution found! 20615673^4=15365639^4+2682440^4+18796760^4  (primitive)
Solution found! 7182177^4=3697823^4+7047520^4+1628600^4
Solution found! 30841113^4=15878887^4+30262880^4+6993400^4
Solution found! 101817921^4=52422079^4+99908960^4+23087800^4
Solution found! 12197457^4=8282543^4+11289040^4+5870000^4  (primitive)
Solution found! 34220961^4=17619039^4+33579360^4+7759800^4
Solution found! 57879897^4=29800103^4+56794720^4+13124600^4
Solution found! 81538833^4=41981167^4+80010080^4+18489400^4
Solution found! 16003017^4=4479031^4+14173720^4+12552200^4  (primitive)
Solution found! 422481^4=217519^4+414560^4+95800^4  (primitive)
Solution found! 16430513^4=16281009^4+3642840^4+7028600^4  (primitive)
Solution found! 47740353^4=24579647^4+46845280^4+10825400^4
Solution found! 37600809^4=19359191^4+36895840^4+8526200^4

Answer 11

I'm not sure how much practical use it is, but as well as d and exactly one of \$a\$, \$b\$, and \$c\$ having to be odd we can exploit modulo 5 and 10 properties.

1. Modulo 5, fourth powers are either 0 or 1. The only way to combine a sum of three of these to make a fourth is with two of the three on the LHS being 0 and with both \$d^4\$ and one of \$a^4\$, \$b^4\$ and \$c^4\$ being 1 modulo 5 (i.e. d ends in 1, 3, 7 or 9 only).

2. In base 10, fourth powers end in either 0000, E1, U6 or 625 where E is an even digit and U is an odd digit. Looking at ways of combining these, and eliminating \$0+0+5=5\$ and \$0+0+6=6\$ as not in their lowest terms, we have that the end digits of the fourth powers satisfy either some permutation of \$0+0+1=1\$ (as in Noam Elkies's solution) or some permutation of \$0+5+6=1\$.

So there must be two multiples of 5 on the LHS and at least one of them must be a multiple of 10 (in Elkies's solution they both are). That should cut down the searching quite a bit too.

I can't see advantage in looking at any other moduli.

There might also be a slight gain in testing merely if the RHS is just a perfect square (a necessary condition of being a fourth power), and only if it meets this test then testing if it is also a square of a square. Any benefit might depend on whether your testing approach to \$d\$ and \$d^4\$ is top down (power to root) or bottom up (root up to its power).

Answer 12

I think I have an \$O(n^2)\$ solution, essentially needing \$O(n^2)\$ space but getting a very large constant factor (~100) shrinkdown. This should make it realistic for the first \$d = 422481\$.

Using Dan Uznanski's odd and even limits

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

(\$a\$ and \$b\$ are even, \$c\$ and \$d\$ are odd)

Using Staycator's mod stuff

\$0+0+1=1\$, \$a^4\$ and \$b^4\$ end in 0000 and \$c\$ and \$d\$ end in 1,3,7 or 9.

Since \$a^4 + b^4 = d^4 - c^4\$, \$c\$ and \$d\$ end in the same four digits.

\$0+6+5=1\$, \$a^4\$ ends in 0000 and \$c^4\$ in 0625. \$b\$ ends in 2, 4, 6 or 8 and d in 1, 3, 7 or 9.

Since \$a^4 + c^4 = d^4 - b^4\$, the last four digits of \$d^4 - b^4\$ are 0625.

Solution

(I haven't put this into code yet, just some notes in Notepad)

• Precalculate \$n^4\$ for \$n\$ 1 to 500k
• Bin by last digit of \$n\$
• Bin \$n^4\$ where n ends in 1, 3, 7 or 9 by last four digits of \$n^4\$
• Bin \$n^4\$ where \$n\$ ends in 2, 4, 6 or 8 by last four digits of \$n^4\$
• Store hash table of \$a^4 + b^4\$ where both \$a\$ and \$b\$ end in zero
• Loop odd last four bins and compare possible \$d^4 - c^4\$ hashes to stored \$a^4 + b^4\$
• Store hash table of \$a^4 + c^4\$ where \$a^4\$ ends in 0000 and \$c^4\$ in 0625
• Find bins where odd \$n\$ minus even \$n\$ bins = 0625 (wrap around). Compare possible \$d^4 - b^4\$ hashes to stored \$a^4 + c^4\$.

Answer 13

Take the idea of calculating \$a^4+b^4\$ and \$d^4-c^4\$ and finding values that are equal. Assume \$a \ge b\$ and \$d \gt c\$.

Create two priority queues, one for the values \$a^4+b^4\$ and one for the values \$d^4-c^4\$, but also keeping track which \$b\$ and \$c\$ were used. Fill the priority queues with values where \$b = 0\$, \$c = d-1\$ and for each \$a, d\$ keep track of which number you have put into the queues (initially \$b = 0\$, \$c = d-1\$). If the smallest number in each queue is the same you have a solution. If not, you remove, for example, (\$a, b\$) from the queue, and insert (\$a, b+1\$) unless \$b+1 \gt a\$; or you remove (\$d, c\$) from the queue and insert (\$d, c-1\$) unless c = 0.

The numbers involved will be big, and 64-bit integers may not be enough. But unless your result get bigger than say \$2^{112}\$ (a, b up to \$2^{28}\$), you can calculate the result once in double precision and once using unsigned 64-bit. The double precision will give you roughly the right result. the unsigned 64 bit will give you the exact last 64 bits of the result. Since \$d, c\$ will be larger, calculate \$d^4 - c^4\$ as \$(d^2 + c^2)(d^2-c^2)\$.

This will take just \$O(n)\$ space, so you'll be unlikely to calculate far enough so this becomes a problem.

Answer 14

In the spirit of actually optimizing the algorithm itself while ignoring the very real overflow problems other have rightly mentioned, think about what is going on here:

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {
            for (d = 1; d < 500000; d++) {
                if (prop(a, b, c, d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

a, b, c, and d are all being raised to the 4th power on every iteration of the innermost loop even though only d has changed.

Also, while addition is a less expensive operation, it adds up (no pun intended) when you do it millions of times for no reason so there's no need to sum the left side of the equation on every iteration either.

This would significantly cut down on calculation time:

for (a = 1; a < 100000; a++) {
    long int a4 = a * a * a * a;
    for (b = 1; b < 300000; b++) {
        long int a4b4 = a4 + b * b * b * b;
        for (c = 1; c < 500000; c++) {
            long int a4b4c4 = a4b4 + c * c * c * c;
            for (d = 1; d < 500000; d++) {
                if (a4b4c4 == d * d * d * d)
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

I haven't taken the time to benchmark this so I can't tell you precisely what kind of improvements you would see but I'll leave that as an exercise for the reader.

Answer 15

With regards to the checking of whether or not \$D\$ is a fourth power, could you not compute, say, all integer 4th powers of \$D\$ up to a limit, and store them in a hash table? That would change the checking time of any potential \$D^4\$ to \$O(1)\$ at an increased prefactor cost at the start. Obviously, you couldn't implement this for an arbitrarily large \$D\$, but OP isn't checking for that (hardcoded limits for \$A,B,C\$), but a ~0.5M sized hashtable is well within the realms of reason.

And if you're checking for bounded \$A,B,C\$, the upper limit for \$D\$ can be found from \$\lceil{\sqrt[4]{A^4+B^4+C^4}}\rceil\$, which would only involve one 4th root operation for the search (validation excluded).

Answer 16

what strategies could I take to make the above code run faster?

Time savers

1. Easy: recognize, at most, combinations, a = 1 to 95800, b = a to 217519, c = b to 414560 and d = c+1 or more need to be tried. Using combinations b = 1 to a-1 and c = 1 to b-1, are redundant. d = 1 to c cannot solve the equation. @phoog

2. After the test values for a & b are set, the values of c and d can be tested very efficiently by incrementing c or d by noticing the lower of sum_ab + c**4 or d**4.

3. Note that if a, b, c, d is a solution, n*a, n*b, n*c, n*d is also a solution. This implies code does not need to check for a solution where a,b,c,d are all even as a solution with an odd value would have been found with smaller values.

4. Every power of 4 of an integer has least significant bits of 0000 or 0001. Adding any 3 of these must equal 0 or 1 (the least significant nibble of d). This implies either all are even (which we can ignore from above) OR 1 of a,b,c must odd, the other 2 even and d must be odd. Using this short-cut reduces the run-time linearly.

5. The task is using integer math exceeding 64 bits. C does not specify a standard type wider than 64-bits. Various platforms have wider types such as __uint128 on gcc.

O(n³) time solution. O(n¹) memory.
The below code exhibits a clear O(n³) trend in time as more and more of the range of a,b,c,d are attempted.

Estimate to complete on my machine a 100% pass: 5,900,000 seconds (69 days).

// Ref 
a = 95800; b = 217519; c = 414560; d = 422481


#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define a (95800 + 1)
#define b (217519 + 1)
#define c (414560 + 1)
#define d (422481 + 1)

// Adjust to some 128 bit type on your platform
// Can leave it at 64-bit if `num/den < 65535/d` or about 0.155 for initial testing
typedef uintmax_t uint128t;

#define pow2(x) ((uint128t)(x)*(x))
#define pow4(x) (pow2((uint_fast64_t)(x)*(x)))
void doit_times(unsigned num, unsigned den) {
  time_t t0;
  time(&t0);
  printf("Portion:%6g%%\t", 100.0 * num / den);
  uint_fast32_t a_max = (uint_fast32_t) round(1.0 * a * num / den);
  uint_fast32_t b_max = (uint_fast32_t) round(1.0 * b * num / den);
  uint_fast32_t c_max = (uint_fast32_t) round(1.0 * c * num / den);
  uint_fast32_t d_max = (uint_fast32_t) round(1.0 * d * num / den);
  printf("a_max:%6" PRIuFAST32 " b_max:%6" PRIuFAST32 //
      " c_max:%6" PRIuFAST32 " d_max:%6" PRIuFAST32 "\t", //
      a_max, b_max, c_max, d_max);
  uint_fast32_t q0 = 1;
  for (uint_fast32_t qa = q0; qa <= a_max; qa++) {
    assert(qa <= 0xFFFF || sizeof(uint128t) > 4);
    uint128t qa4 = pow4(qa);
    for (uint_fast32_t qb = qa; qb <= b_max; qb++) {
      uint128t sum_ab4 = qa4 + pow4(qb);
      uint_fast32_t qc = qb;
      uint128t sum_abc4 = sum_ab4 + pow4(qc);
      uint_fast32_t qd = qc + 1;
      uint128t qd4 = pow4(qd);

      // Loop until `c` to too large.
      for (;;) {
        if (sum_abc4 < qd4) {
          if (qc >= c_max) break;
          qc++;
          sum_abc4 = sum_ab4 + pow4(qc);
        } else if (sum_abc4 > qd4) {
          //if (qd >= d_max) break;
          qd++;
          qd4 = pow4(qd);
        } else {
          puts("Success!");
          printf("a %" PRIuFAST32 "\n", qa);
          printf("b %" PRIuFAST32 "\n", qb);
          printf("c %" PRIuFAST32 "\n", qc);
          printf("d %" PRIuFAST32 "\n", qd);
          fflush(stdout);
          exit(0);
        }
      }

    }
  }
  time_t t1;
  time(&t1);
  double dt = difftime(t1, t0);
  printf("t:%g\n", dt);
  fflush(stdout);
}

int main(void) {
  for (unsigned num = 1; num <= 1000; num++) {
    doit_times(num, 1000);
  }
  puts("Done");
  return -1;
}

Output

Portion:   0.1% a_max:    96 b_max:   218 c_max:   415 d_max:   422 t:0
Portion:   0.2% a_max:   192 b_max:   435 c_max:   829 d_max:   845 t:0
Portion:   0.3% a_max:   287 b_max:   653 c_max:  1244 d_max:  1267 t:1
Portion:   0.4% a_max:   383 b_max:   870 c_max:  1658 d_max:  1690 t:1
Portion:   0.5% a_max:   479 b_max:  1088 c_max:  2073 d_max:  2112 t:2
Portion:   0.6% a_max:   575 b_max:  1305 c_max:  2487 d_max:  2535 t:3
Portion:   0.7% a_max:   671 b_max:  1523 c_max:  2902 d_max:  2957 t:5
Portion:   0.8% a_max:   766 b_max:  1740 c_max:  3316 d_max:  3380 t:6
Portion:   0.9% a_max:   862 b_max:  1958 c_max:  3731 d_max:  3802 t:10
Portion:     1% a_max:   958 b_max:  2175 c_max:  4146 d_max:  4225 t:12
Portion:   1.1% a_max:  1054 b_max:  2393 c_max:  4560 d_max:  4647 t:16
Portion:   1.2% a_max:  1150 b_max:  2610 c_max:  4975 d_max:  5070 t:21
Portion:   1.3% a_max:  1245 b_max:  2828 c_max:  5389 d_max:  5492 t:25
Portion:   1.4% a_max:  1341 b_max:  3045 c_max:  5804 d_max:  5915 t:32
Portion:   1.5% a_max:  1437 b_max:  3263 c_max:  6218 d_max:  6337 t:39
Portion:   1.6% a_max:  1533 b_max:  3480 c_max:  6633 d_max:  6760 t:48
Portion:   1.7% a_max:  1629 b_max:  3698 c_max:  7048 d_max:  7182 t:56
Portion:   1.8% a_max:  1724 b_max:  3915 c_max:  7462 d_max:  7605 t:67
Portion:   1.9% a_max:  1820 b_max:  4133 c_max:  7877 d_max:  8027 t:77
Portion:     2% a_max:  1916 b_max:  4350 c_max:  8291 d_max:  8450 t:91
Portion:   2.1% a_max:  2012 b_max:  4568 c_max:  8706 d_max:  8872 t:104
Portion:   2.2% a_max:  2108 b_max:  4785 c_max:  9120 d_max:  9295 t:120

Note that testing d is within a limit is not needed as c exceeds the limit first or shortly thereafter.

Answer 17

How about finding the powers first and then computing the equivalency? This will eliminate the recomputation of the fourth power of the same number (due to space considerations).

Something like this can be programmed:

%from numbers 1 to 500000
for k=1:50000:500000
%for indexing
i=0;
%compute the forth powers from k to k+49999 and store the powers in a 
for j=k:k+50000-1
a[i++]=j*j*j*j;
end
for x=1:50000
for y=1:50000
for z=1:50000
for t=1:50000
  if (a[x]+a[y]+a[z]==a[t])
   return;
end
end
end
end

Other optimization methods can also be used.