135
\$\begingroup\$

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?

\$\endgroup\$
  • 16
    \$\begingroup\$ one thing you can do is for (b = a; ... for (c = b;. For example, you don't need to check both (a, b, c) = (1, 2, 3) and (a, b, c) = (3, 2, 1). Similarly, for any given set of three values for (a, b, c), it's extremely wasteful to check all values of d between 1 and 500,000. Instead, calculate d^4 by adding the fourth powers of (a, b, c) and then check whether its fourth root is an integer. This will improve the absolute time necessary, but not the general time complexity. \$\endgroup\$ – phoog Oct 25 '16 at 18:11
  • 12
    \$\begingroup\$ You can learn a lot by studying the improvements suggested in the answers - but none of those answers turns your impractical method into one that you could actually use to find Elkies's example. You need much more than good programming. That's why it took a mathematician like Elkies to come up with it. \$\endgroup\$ – Ethan Bolker Oct 26 '16 at 13:26
  • 2
    \$\begingroup\$ It is the case (~1e20 ops) when a better algorithm (better Big O—time complexity) is more important than simple implementation tricks that improve the implementation by a constant factor (even 100 times faster is too slow in this case). \$\endgroup\$ – jfs Oct 26 '16 at 18:02
  • 2
    \$\begingroup\$ There are 100000 * 300000 * 500000, or 1.5e16 values within your ranges for a b and c. Even if you could check whether a number is a perfect 4th power in 1 nanosecond, it would take 174 years to explore the entire search space. So the most important thing is going to be drastically cutting down on the number of candidates to search somehow. \$\endgroup\$ – Tavian Barnes Oct 26 '16 at 20:10
  • 4
    \$\begingroup\$ @diego.martinez: pre or post increment are completely equivalent here, any decent compiler will generate the same code, including if one wrote a = a + 1 instead of a++ \$\endgroup\$ – chqrlie Oct 28 '16 at 21:27

17 Answers 17

136
\$\begingroup\$

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.

\$\endgroup\$
  • 8
    \$\begingroup\$ Integral pow4 is perfectly reasonable, but you need to use a datatype large enough to hold the result! \$\endgroup\$ – Ben Voigt Oct 25 '16 at 20:26
  • 7
    \$\begingroup\$ There's one major problem with switching to integer arithmetic: the solution involves the number 31858749840007945920321, which cannot be stored in even a 64-bit int. \$\endgroup\$ – Mark Oct 25 '16 at 21:23
  • 55
    \$\begingroup\$ @Mark: And therefore it cannot be stored in a 64 bit double either. \$\endgroup\$ – Eric Lippert Oct 25 '16 at 22:45
  • 3
    \$\begingroup\$ @ymbirtt int multiplication is properly associative (ignoring overflow as the compiler is allowed to do for signed) so the compiler can do that optimization. \$\endgroup\$ – ratchet freak Oct 26 '16 at 8:26
  • 10
    \$\begingroup\$ @ymbirtt: GCC will optimize x*x*x*x into two imul instructions: pow4: movq %rdi, %rax; imulq %rdi, %rax; imulq %rax, %rax; ret (those will be inlined later again, by the way). So yes, the compiler is intelligent enough for this optimization :). \$\endgroup\$ – Zeta Oct 26 '16 at 8:34
76
\$\begingroup\$

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 d4 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 a4 + b4 + c4 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 1284 equal to sum? No. Is it bigger? Yes. Then next try 644. Is it equal to sum? No. Is it smaller? Yes. Then try 964, 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.

\$\endgroup\$
  • \$\begingroup\$ Which would be faster? A binary search or the sqrt(sqrt(pow(pow(x, 2), 2))) that you suggest first? I suppose the answer is for me to implement both and see. Great answer. Thanks man. \$\endgroup\$ – Aidenhjj Oct 25 '16 at 19:14
  • 7
    \$\begingroup\$ @Aidenhjj: Eric isn't suggesting sqrt(sqrt(pow(pow(x, 2), 2))) but pow(pow(sqrt(sqrt(x)), 2), 2). The order is extremely important. \$\endgroup\$ – Ben Voigt Oct 25 '16 at 20:22
  • 2
    \$\begingroup\$ A sqrt implementation likely uses binary search internally so doing two sqrts instead of one binary search probably won't be a win. \$\endgroup\$ – 5gon12eder Oct 26 '16 at 0:12
  • 2
    \$\begingroup\$ Rather than binary search for a d whose fourth power is equal to the sum of a⁴+b⁴+c⁴ in general, it might be worth remember the value of d that was just larger than the sum (call it dlast), then the next time round the loop for c check if this is now equal to the sum (win!), still greater (leave it alone), or now less than the sum (in which case increment by one and try again). That should require only one or two tries in the inner loop (rather than O(log(n))) \$\endgroup\$ – Martin Bonner Oct 26 '16 at 10:35
  • 5
    \$\begingroup\$ Another comment: It's probably worth storing a4 == pow4(a) and a4b4 == a4 + pow4(b). Still need to make O(n³) calls to pow4 for c and d but it still saves approximately half the calls. \$\endgroup\$ – Martin Bonner Oct 26 '16 at 10:38
31
\$\begingroup\$

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.

\$\endgroup\$
  • 5
    \$\begingroup\$ If you're running on a 64-bit system, GCC supports an __int128 type that can handle the numbers needed. \$\endgroup\$ – Mark Oct 25 '16 at 21:27
  • 2
    \$\begingroup\$ You don't need arbitrary-precision, just double-word math, which is nowhere near as slow. Keeping just the low bits as a fast filter is a good idea, but it breaks monotonicity so the binary-search method won't work anymore, neither will any other sqrt function. But it should apply beautifully to @vnp's sieve-based search. \$\endgroup\$ – Ben Voigt Oct 25 '16 at 22:13
  • 3
    \$\begingroup\$ A long int in C has at least 32 bits. \$\endgroup\$ – 5gon12eder Oct 26 '16 at 0:14
30
\$\begingroup\$

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.

\$\endgroup\$
  • 17
    \$\begingroup\$ Actually it gets even better than this: all even numbers, when taken the to the fourth power, are divisible by 16, and all odd numbers, when taken to the fourth power, have a remainder of 1 when divided by 16. So it must be -- since having all of a, b, c, and d even means there's a smaller solution -- that d must be odd and exactly one of a, b, and c must be odd as well, because that's the only way to get the same remainder mod 16 on both sides. \$\endgroup\$ – Dan Uznanski Oct 26 '16 at 1:42
  • \$\begingroup\$ So with @DanUznanski simplification, if a is odd, both b and c must be even, if a is even, b can be either, and c must be the opposite parity to b.. That means the search space is cut in just over half. Can we do similar tricks with mod 3 calculations? \$\endgroup\$ – Martin Bonner Oct 26 '16 at 10:52
  • 1
    \$\begingroup\$ "since if they are all even the gcd is 2"? - 2 isn't the gcd. \$\endgroup\$ – JimmyB Oct 27 '16 at 14:02
  • \$\begingroup\$ @JimmyB true, not necessarily. I'll edit to say common divisor. \$\endgroup\$ – mascoj Oct 27 '16 at 14:56
  • 3
    \$\begingroup\$ @MartinBonner -- Taking modulo 3 says that either none of a, b, and c can be multiples of 3, or that exactly two of a, b, and c must be multiples of 3. Modulo 5 is a bit more helpful: Exactly two of a, b, and c must be multiples of 5. \$\endgroup\$ – David Hammen Oct 31 '16 at 15:55
30
\$\begingroup\$

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);
    }
  }

}
\$\endgroup\$
  • \$\begingroup\$ Can you explain to me the syntax of auto sumcmp =[&list,&sumptr](int a,int b){ return list[sumptr[a]]+list[a]>list[sumptr[b]]+list[b]; };. Is that a function definition or are you defining a new vector? \$\endgroup\$ – Aidenhjj Oct 27 '16 at 9:29
  • \$\begingroup\$ I guess you're defining a new vector where you loop through list in var a and loop through simultaneously sumptr in var b. But then I get a bit lost. Here, do a and b have the same meaning as in my question? Because the way they're treated here, from my limited understanding it seems that a is really a^4 and b is b. \$\endgroup\$ – Aidenhjj Oct 27 '16 at 9:44
  • 1
    \$\begingroup\$ Can you make your naming conventions a little more explicit maybe? \$\endgroup\$ – Aidenhjj Oct 27 '16 at 9:45
  • \$\begingroup\$ Aidenhjj: I've changed the variable names to match the original problem description. auto sumcmp = (etc) is a so-called lambda function definition which behaves like a function but can reference the variables list and sumptr. Sorry, this is probably a distraction from the main purpose. It could be avoided using global variables and functions, or a slightly different data structure. \$\endgroup\$ – Alex Selby Oct 27 '16 at 13:33
  • 1
    \$\begingroup\$ @AlexSelby "you need to find the value of \$a\$ such that \$a^4 + b_a^4\$ is the smallest" I don't understand what you mean here. Isn't that just \$a = 1\$ always? \$\endgroup\$ – Tavian Barnes Oct 27 '16 at 15:48
21
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Good point about Frye, but you don't need any fancy stuff you mentioned to get it down to hours on a modern machine: simple brute force yields the answer in a few hours. If you put some effort into making it fast with modulo restrictions, hand-coded (rather than library) data structures, multithreading, etc., then it would probably run in a few minutes or less. \$\endgroup\$ – Alex Selby Oct 30 '16 at 11:01
  • \$\begingroup\$ @Edward: I posted an example (unoptimised) program with my answer which finds a solution in a few hours. \$\endgroup\$ – Alex Selby Oct 31 '16 at 23:10
  • \$\begingroup\$ I posted a different method with some optimisations, which finds the solution in about 40 seconds. \$\endgroup\$ – Alex Selby Nov 2 '16 at 1:11
18
\$\begingroup\$

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})\$.

\$\endgroup\$
  • 7
    \$\begingroup\$ Building a table of sums takes \$O(n^2)\$ time AND SPACE, which may limit its applicability. \$\endgroup\$ – Ben Voigt Oct 25 '16 at 20:25
  • 2
    \$\begingroup\$ @BenVoigt Correct; however there is a way with an \$O(n)\$ space and \$O(n^2)\$ time. I intentionally didn't get into it. \$\endgroup\$ – vnp Oct 25 '16 at 21:04
  • 2
    \$\begingroup\$ @JS1 At the afterthought I am curious myself. Give me by the end of the day, and either I post a code, or admit myself wrong. \$\endgroup\$ – vnp Oct 25 '16 at 21:48
  • 9
    \$\begingroup\$ @JS1 Day is over. I didn't provide the code, so I must admit myself wrong. \$\endgroup\$ – vnp Oct 26 '16 at 7:43
  • 1
    \$\begingroup\$ @gnasher729: But no one cares about n^4 operations, since several answers have shown algorithms in O(n^3) time and O(1) space. \$\endgroup\$ – Ben Voigt Oct 27 '16 at 19:04
12
\$\begingroup\$

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);
}
\$\endgroup\$
  • 1
    \$\begingroup\$ Hate the formatting, but the code is solid. You can save even more time by noting that a and b can't both be odd. If you add add one more check before list0.push_back(...), like this: if (!(a & 1 & b)) then the threaded version saves about 25% of the time on my machine (69 seconds vs. 90 for the original). \$\endgroup\$ – Edward Nov 2 '16 at 2:55
  • \$\begingroup\$ Good point. I wasn't going to get into other moduli improvements as they get quite complicated, but your improvement is simple and certainly worth having. I've edited the code above to incorporate the change. (Sorry about the formatting - it's what I'm used to.) \$\endgroup\$ – Alex Selby Nov 2 '16 at 9:43
  • \$\begingroup\$ I thought about accepting this as the answer, but Zeta's answer is very complete and has helped me learn a lot. As and when I've picked up enough to fully understand your answers, I'll come back to them. \$\endgroup\$ – Aidenhjj Nov 2 '16 at 12:39
  • \$\begingroup\$ I think it was my little joke that made it popular! Otherwise, it's a pretty standard 'I don't know anything about programming' question. \$\endgroup\$ – Aidenhjj Nov 2 '16 at 12:59
  • \$\begingroup\$ @Aidenhjj No problem - I'm not worried about this answer being accepted or not. I think it is legitimate to interpret your question as "explain what programming techniques would improve this specific method?", which Zeta did, or "how might you do this better if you can change the algorithm a bit?" which I chose to do. Your question brings up quite a lot of useful techniques that occur a lot. The meet-in-the-middle technique is applicable to many search situations, reducing a search of \$N\$ to \$\sqrt{N}\$. Binning into independent sets is often used to reduce memory requirements. Etc. \$\endgroup\$ – Alex Selby Nov 2 '16 at 18:09
8
\$\begingroup\$

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();
}
\$\endgroup\$
  • \$\begingroup\$ Estimated run-time: 5 yrs, \$\endgroup\$ – chux Oct 26 '16 at 2:42
  • \$\begingroup\$ It's not entirely clear what the code is actually doing mathematically. Could you explain it? \$\endgroup\$ – Edward Oct 26 '16 at 3:17
  • \$\begingroup\$ @Edward What parts of code were not clear? \$\endgroup\$ – chux Oct 26 '16 at 15:20
  • \$\begingroup\$ uintmax_t is not 128 bits on any platform I know of, certainly not x86_64 Linux with gcc. You'll need to use __uint128_t. \$\endgroup\$ – Tavian Barnes Oct 26 '16 at 15:26
  • \$\begingroup\$ @TavianBarnes The // Adjust to 128 bit type typedef uintmax_t uint_fast128t; is simply a place-holder for whatever 128-bit type a compiler may have access - adjust as needed. __uint128_t is not part of standard C, yet an efficient solution would use some extension such as __uint128_t . IAC, larger than 64-bit math is not needed until factor < 7. \$\endgroup\$ – chux Oct 26 '16 at 15:30
7
\$\begingroup\$

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).

\$\endgroup\$
7
\$\begingroup\$

The solution you are talking about was not found by Elkies but by Fry. 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 theorem 1 and two can be found here and a summary that is referenced from here 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.

At the moment (2016-11-04) the file erroneously contains the program twice. Delete all lines from line 1478 to the end.

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
here is a summary and here is the first page of Wards's article


[Robert Gerbicz' euler413.c] https://sites.google.com/site/robertgerbicz/euler413.c

[Robert Gerbicz' euler413.exe] https://sites.google.com/site/robertgerbicz/euler413.exe



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 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
\$\endgroup\$
  • \$\begingroup\$ Unfortunately, the cited euler413.c program is malformed and does not compile. \$\endgroup\$ – Edward Nov 4 '16 at 15:39
  • \$\begingroup\$ @Edward: That is true. The file erroneously contains the program twice. Delete everything from line 1438 until the last line. \$\endgroup\$ – miracle173 Nov 4 '16 at 16:45
  • \$\begingroup\$ @Edward there is also a compiled version euler413.exe for windows on Robert Gerbicz' page. I added some notes in the text. \$\endgroup\$ – miracle173 Nov 4 '16 at 17:49
  • \$\begingroup\$ I appreciate the thought, but line 1438 is the last line in what I downloaded and I have a policy of not running random executables from the internet for safety reasons. \$\endgroup\$ – Edward Nov 4 '16 at 18:29
  • \$\begingroup\$ The code compiles for me, but to make it 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). \$\endgroup\$ – Alex Selby Nov 4 '16 at 23:55
6
\$\begingroup\$

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\$.
\$\endgroup\$
4
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Initially you have the smallest possible values of \$a^4 + b^4\$ but the largest possible values of \$d^4 - c^4\$ -- you need some serious refinement to this approach before it will let you start eliminating any pairs from consideration. \$\endgroup\$ – Ben Voigt Oct 27 '16 at 19:08
3
\$\begingroup\$

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 on 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.

\$\endgroup\$
2
\$\begingroup\$

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(n3) time solution. O(n1) memory.
The below code exhibits a clear O(n3) 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.

\$\endgroup\$
0
\$\begingroup\$

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).

\$\endgroup\$
  • \$\begingroup\$ If the loops are ordered such that a >= b >= c, d is in the range [a ; pow(3, 0.25) * a]. \$\endgroup\$ – chqrlie Oct 28 '16 at 22:33
  • \$\begingroup\$ Good point. Saves a bit of memory. \$\endgroup\$ – M. Withnall Oct 29 '16 at 0:33
-1
\$\begingroup\$

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.

\$\endgroup\$
  • 1
    \$\begingroup\$ It doesn't seem to improve neither a memory usage nor an expected run time of the original code. \$\endgroup\$ – CiaPan Oct 28 '16 at 7:35
  • \$\begingroup\$ Why? In the original code the author calculates the fourth power of the same numbers four times, .e.g. for a, b, c and d, separately. The above code calculates the fourth power just once and use the calculated values when necessary. There may not be a memory optimization. But my priority is speed rather than the memory here. \$\endgroup\$ – tempx Nov 1 '16 at 6:15
  • 1
    \$\begingroup\$ Your pre-computing of 4-th powers doesnt save you much. Your code is going to perform \$5000^4 = 625\cdot 10^{16}\$ iterations, just like the original code. Even if you make each iteration time \$1/100\$ of the original, the whole execution will still take longer than anyone could wait. See the comment by Tavian Barnes. \$\endgroup\$ – CiaPan Nov 1 '16 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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