Disproving Euler proposition by brute force in C 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 &lt;stdio.h&gt;
#include &lt;time.h&gt;
#include &lt;math.h&gt;

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 &lt; 100000; a++) {
        for (b = 1; b &lt; 300000; b++) {
            for (c = 1; c &lt; 500000; c++) {
                for (d = 1; d &lt; 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?
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&shy;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 &lt;stdio.h&gt;
#include &lt;time.h&gt;
#include &lt;math.h&gt;


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 &lt; 100000; a++) 
    for (b = 1; b &lt; 300000; b++) 
        for (c = 1; c &lt; 500000; c++) 
            for (d = 1; d &lt; 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 &lt; 100000; a++) 
    for (long int b = 1; b &lt; 300000; b++) 
        for (long int c = 1; c &lt; 500000; c++) 
            for (long int d = 1; d &lt; 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&shy;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 &lt; 100000; a++) {
        a4 = pow4(a);                             // remember
        for (b = a; b &lt; 300000; b++) {
            b4 = pow4(b);                         // remember
            for (c = b; c &lt; 500000; c++) {
                c4 = pow4(c);                     // the fourth power
                for (d = c; d &lt; 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&shy;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 &lt; 100000; a++) {
        const Number a4 = pow4(a);                         // remember
        for (Number b = a; b &lt; 300000; b++) {
            const Number b4 = pow4(b);                     // remember
            for (Number c = b; c &lt; 500000; c++) {
                const Number c4 = pow4(c);                 // the fourth power
                for (Number d = c; d &lt; 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.