This problem was posed to an acquaintance of mine (at their university) and piqued my interest so I tried to solve it. The description goes as follows:
Write a program that finds a solution to the following equation: $x^4=a^4+b^4+c^4+d^4$. One solution is $651^4 = 240^4 + 340^4 + 430^4 + 599^4$ however you are tasked to finding a smaller solution. The program is not allowed to use any type of array/matrix.
The goal is to solve this in a reasonable/fastest amount of time. I have designed a solution that finishes in 10 seconds starting from $x=2$ up to $x=353$ (I know runtimes are arbitrary per machine but I don't want to do an asymptotic analysis at the moment). However I'm not at all versed in this kind of mathematics so I'm looking to see if there's a different approach or some algorithmic optimization I'm not aware of that would make this much faster.
A few things to note:
- The weird array/matrix restriction essentially blocks dynamic programming. Not sure why they added that but w/e.
- The smaller solution is $353^4 = 30^4 + 120^4 + 272^4 + 315^4$. This might will help experimentation.
- The specific problem was to be solved in C language but here I'm looking for an algorithm so pseudo code or any language (no esoteric languages please, not even stuff like APL) is fine.
- No multithreading/parallelism.
- A good solution should be able to find the 3rd biggest solution in a reasonable amount of time $2487^4 = 435^4 + 710^4 + 1384^4 + 2420^4$.
EDIT: Was asked to describe my current solution. The basic premise is a brute force one where I run incrementing loops for each value of $x$ followed by $a, b, c$ etc. There are two "tricks".
- Each successive variable is always larger than the previous one $d > c > b > a$ and all are smaller than $x-1$.
- Since $d$ is the largest and thus would take the longest to iterate through, I don't iterate but simply check whether $\sqrt[4]{x^4 - a^4 - b^4 - c^4} = d$ is an integer.
Here is some fully functional C code:
#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>
/**
* Checks if the provided number n is a perfect square before
* taking the actual square root (which is a very slow operation).
*/
uint64_t perfectSquareRoot(uint64_t n)
{
uint64_t h = n & 0xF; // last hexadecimal "digit"
if (h > 9)
return 0; // return immediately in 6 cases out of 16.
// Take advantage of Boolean short-circuit evaluation
if ( h != 2 && h != 3 && h != 5 && h != 6 && h != 7 && h != 8 )
{
uint64_t sqrtN = sqrt(n);
uint64_t t = (uint64_t) floor(sqrtN + 0.5 );
if(t*t == n) {
return sqrtN;
}
}
return 0;
}
int main(int argc, char** argv) {
int x,a,b,c,d;
clock_t begin = clock();
for(x = 2; x < 651; x++) {
printf("Trying x = %d\n", x);
uint64_t x4 = pow(x,4);
for(a=1; a<x-1; a++) {
uint64_t a4 = pow(a,4);
uint64_t ab4 = 0;
for(b=a; ab4 <= x4; b++) {
ab4 = a4+pow(b,4);
uint64_t abc4 = 0;
for(c=b; abc4 <= x4; c++) {
abc4 = ab4+pow(c,4);
uint64_t d4 = x4 - abc4; // take the difference between x^4 - a^4+b^4+c^4 = d^4
uint64_t d2 = perfectSquareRoot(d4); // if d4 is a perfect square of a perfect square we found our number
if(d2) {
d = perfectSquareRoot(d2);
if(d) {
if(x4 == abc4 + pow(d,4)) { // do a final reverse check to handle numerical instability
printf("x=%d, a=%d, b=%d, c=%d, d=%d\n", x,a,b,c,d);
goto endlbl;
}
}
}
}
}
}
}
endlbl:;
clock_t end = clock();
printf("Time Spent: %f", (double) (end - begin / CLOCKS_PER_SEC));
return 0;
}
Part one, brute force: Here is an alternative brute force implementation. The code is julia code, copy+pasted in Appendix 1, my first attempt to implement something in julia, i took this approach after having the "same code" in sage and being unhappy with the waiting time. Since julia should have a quick way to go with loops & CO, i gave it a try.
Some words on the implementation. Consider solutions of the equation: $$ \tag{$*$} x^4 = a^4+b^4+c^4+d^4\ , $$ and try to enumerate them. (This is different from finding a(n infinite) family of solutions of the above diophantine equation, instead we want to list all possible $(x;a,b,c,d)$ below some bound, and check that we have $(*)$. In the process, we can of course eliminate "in time" some values, if there is a good reason to do this.) By symmetry, we can restrict to those variables that satisfy: $$ x > a \ge b\ge c\ge d>0\ . $$ In particular, we have: $$ \begin{aligned} a^4 &<x^4\underbrace{\le}_{\text{but }\ne} a^4+a^4+a^4+a^4=4a^4\ ,&&\text{so} & \frac 1{4^{1/4}}x &<a<x\ , \\ N &:= x^4 -a^4=:n^4\ ,\\ n^4 = N & = b^4 + c^4 + d^4 < b^4+b^4+b^4=3b^4\ ,&&\text{so} & \frac 1{3^{1/4}}n &\le b\le\min(a, n)\ , \\ M &:= x^4 -a^4 -b^4=:m^4\ ,\\ m^4 = M & = c^4 + d^4 < c^4+c^4=2c^4\ ,&&\text{so} & \frac 1{2^{1/4}}m &\le c\le\min(b, m)\ . \end{aligned} $$ With these restrictions the hierarchic loops for
x, thena, thenb, thencare thiner. The brute force does not spent time at irrelevant places. We get the following results in time: $$ \begin{array}{|r||r|r|r|r|l|} \hline x&a&b&c&d&\text{time in between}\\\hline\hline 353 & 315 & 272 & 120 & 30 & \text{after }2.819s\\\hline 651 & 599 & 430 & 340 & 240 & \text{after further }17.146s\\\hline 2487 & 2420 & 1384 & 710 & 435 & \text{after further }4200.359s\sim1h\; 10'\;1'' \\\hline 2501 & 2365 & 1432 & 1190 & 1130 & \text{after further }90.856s\sim1'\; 30'' \\\hline 2829 & 2745 & 1546 & 1010 & 850 & \text{after further }2837.045s\sim47'\; 17'' \\\hline 3723 & 3152 & 2460 & 2345 & 2270 & \text{after further }15891.990s\sim4h\;24'\; 51'' \\\hline 3973 & 3395 & 3230 & 1652 & 350 & \text{after further }5967.016s\sim1h\;39'\; 27'' \\\hline 4267 & 4094 & 2650 & 1060 & 205 & \text{after further }9992.362s\sim2h\;45'\; 32'' \\\hline \end{array} $$ This was the brute force approach, using it i've got the best results. Minimal brain, maximal success. One can already stop here. (But since the question asks also for ideas to proceed, here is some added information about what i also tried...)Part two, using the arithmetic of $\Bbb Z[i]$:
I have also tried to restrict the part of the search by means of arithmetics. To illustrate the situation, i will pick some random value for $(x,a,b)$, and show how to search for $(c,d)$ in a different manner, not looping through all possible values. First of all, we take some random case where this may make sense:
(A $c$-loop is now maybe not a too good idea.) The question, the partial programming task is as follows, we have some numbers $x,a,b$. Can we find in this special case some $c,d$ satisfying $(*)$?
This means: $$c^4+d^4= M:=32961715639129459861921278445069934729992837454625353749504\ .$$ Let $C=c^2$, $D=d^2$. Can we find at least $C,D$ first, all possible choices?! Yes, and we can do this algorithmically as follows. If $C^2+D^2=M$, then we have in the gaussian ring $\Bbb Z[i]$ (with unique factorization): $$ M =(C+iD)(C-iD)\ . $$ Here $i=\sqrt{-1}$. Each factorization over $\Bbb Z[i]$ comes from one over $\Bbb Z$, so the first step is to factor $M$:
And since the factor $7$, a prime congruent to three modulo four, comes to an odd power, we can already eliminate this case from the discussion, there is no chance to factorize. We have a quick verdict after a long/expensive factorization process.
OK, we are generating so long random cases, till we get only even powers for such "bad" (inert) primes.
We have a big candidate $M$, and try to write it in the form $C^2+D^2$ first. All ways to do this are needed. For this, we write the prime pieces in this factorization as wanted, we factorize them over $\Bbb Z[i]$, in code $i$ is
jto not distroy the default $i$ variable in sage:(Here, $\pm 1$, $\pm i$ are units, we ignore them for our purpose. Also $1\pm i$ are conjugated, we will place $1+i$ on one side, $1-i$ on the other side, but it doesn't matter, they are conjugated...) And now we separate in conjugated pieces of same norm, each case with each case, and build products. All the needed products are, taking $(1+2i)$ on a fixed side: $$ \small \begin{aligned} &(1+2i)(1+i)^4(39050+10097i)(24613485632209445972735+12375789567144879751684i) &&\text{and}\\ &(1-2i)(1-i)^4(39050-10097i)(24613485632209445972735-12375789567144879751684i)\\ \\[3mm] &(1+2i)(1+i)^4(39050+10097i)(24613485632209445972735-12375789567144879751684i) &&\text{and}\\ &(1-2i)(1-i)^4(39050-10097i)(24613485632209445972735+12375789567144879751684i)\\ \\[3mm] &(1+2i)(1+i)^4(39050-10097i)(24613485632209445972735+12375789567144879751684i) &&\text{and}\\ &(1-2i)(1-i)^4(39050+10097i)(24613485632209445972735-12375789567144879751684i)\\ \\[3mm] &(1+2i)(1+i)^4(39050-10097i)(24613485632209445972735-12375789567144879751684i) &&\text{and}\\ &(1-2i)(1-i)^4(39050+10097i)(24613485632209445972735+12375789567144879751684i)\\ \end{aligned} $$ We multiply, and get for the first case:
Now we would adjust signs (for
Donly) and look if $C=c^2$ and $D=d^2$. At any rate, the factorization was replacing the $c$-loop, which becomes an issue only for "bigger" $M$-values involved.How does it work with a solution? Consider the triple $x,a,b=2829,2745,1546$. Then $$ M=2829^4-2745^4-1546^4= 2^5 \cdot 5^4 \cdot 73 \cdot 241 \cdot 4441 \\ =(1+i)^5(1-i)^5\cdot(2+i)^4(2-i)^4\cdot(8+3i)(8-3i)\cdot(15+4i)(15-4i) \cdot (60+29i)(60-29i)\ . $$ Now we consider from each pair only one number, all the list is, fixing $(60+29i)$ "on one side": $$ \begin{aligned} (1+i)^5\cdot(2+i)^4(2-i)^0\cdot(8+3i)\cdot(15+4i)\cdot(60+29i) &= 672452i + 1053764 \ ,\\ (1+i)^5\cdot(2+i)^3(2-i)^1\cdot(8+3i)\cdot(15+4i)\cdot(60+29i) &= -439540i + 1170220 \ ,\\ (1+i)^5\cdot(2+i)^2(2-i)^2\cdot(8+3i)\cdot(15+4i)\cdot(60+29i) &=-1199900i + 350500 \ ,\\ (1+i)^5\cdot(2+i)^1(2-i)^3\cdot(8+3i)\cdot(15+4i)\cdot(60+29i) &=-1000340i - 749620 \ ,\\ (1+i)^5\cdot(2+i)^0(2-i)^4\cdot(8+3i)\cdot(15+4i)\cdot(60+29i) &= -508i - 1250044 \ ,\\[3mm] (1+i)^5\cdot(2+i)^4(2-i)^0\cdot(8-3i)\cdot(15+4i)\cdot(60+29i) &=-186244i + 1236092 \ ,\\ (1+i)^5\cdot(2+i)^3(2-i)^1\cdot(8-3i)\cdot(15+4i)\cdot(60+29i) &= -1100620i + 592660 \ ,\\ (1+i)^5\cdot(2+i)^2(2-i)^2\cdot(8-3i)\cdot(15+4i)\cdot(60+29i) &= -1134500i - 524900 \ ,\\ (1+i)^5\cdot(2+i)^1(2-i)^3\cdot(8-3i)\cdot(15+4i)\cdot(60+29i) &= -260780i - 1222540 \ ,\\ (1+i)^5\cdot(2+i)^0(2-i)^4\cdot(8-3i)\cdot(15+4i)\cdot(60+29i) &= 821564i - 942148 \ ,\\[3mm] (1+i)^5\cdot(2+i)^4(2-i)^0\cdot(8+3i)\cdot(15-4i)\cdot(60+29i) &= 58468i + 1248676 \ ,\\ (1+i)^5\cdot(2+i)^3(2-i)^1\cdot(8+3i)\cdot(15-4i)\cdot(60+29i) &= -963860i + 795980 \ ,\\ (1+i)^5\cdot(2+i)^2(2-i)^2\cdot(8+3i)\cdot(15-4i)\cdot(60+29i) &=-1215100i - 293500 \ ,\\ (1+i)^5\cdot(2+i)^1(2-i)^3\cdot(8+3i)\cdot(15-4i)\cdot(60+29i) &=-494260i - 1148180 \ ,\\ (1+i)^5\cdot(2+i)^0(2-i)^4\cdot(8+3i)\cdot(15-4i)\cdot(60+29i) &=621988i - 1084316 \ ,\\[3mm] (1+i)^5\cdot(2+i)^4(2-i)^0\cdot(8-3i)\cdot(15-4i)\cdot(60+29i) &=-776996i + 979228 \ ,\\ (1+i)^5\cdot(2+i)^3(2-i)^1\cdot(8-3i)\cdot(15-4i)\cdot(60+29i) &=-1249580i - 34060 \ ,\\ \color{blue}{ (1+i)^5\cdot(2+i)^2(2-i)^2\cdot(8-3i)\cdot(15-4i)\cdot(60+29i) } & \color{blue}{=-722500i - 1020100} \ ,\\ (1+i)^5\cdot(2+i)^1(2-i)^3\cdot(8-3i)\cdot(15-4i)\cdot(60+29i) &=382580i - 1190060 \ ,\\ (1+i)^5\cdot(2+i)^0(2-i)^4\cdot(8-3i)\cdot(15-4i)\cdot(60+29i) &=1181596i - 407972 \ , \end{aligned} $$ In each of the cases we have to adjust the sign of the real, and of the imaginary part, then check if both are squares of integers. The solution comes from the blue line. If there are any "inert" primes (those which are three modulo four) in the decomposition, they have to be distributed equally on each side.
This explicit procedure may be also used to cover some cases, it may be useful when we search for solutions for "higher" $x$ values in some small range.
Moral: The factorization part is the weak step, so this solution is more qweak then quick.
Appendix 1: Used Julia code in the sequel. Some comments in advance. The package
Heckeis imported only to have thatis_squarefunctionality. There is no point for me in making this function even more performant. I am happy to wait those seconds or hours when the code runs, my focus is on a (mathematically) structural way to implement the objects. The clock is given byTickTock, when typingtick()it starts ticking, at atock()it stops, delivering info to the amount of time passed. So below, we first load some libraries, this takes some good seconds, after that the clock ticks, and in some few seconds we have the first solution. We stop the clock, show the solution, check it, then start again the clock from zero, the second solution comes then in less than half a minute after the first one. Again we stop the clock, show the solution, check it, and start hunting the third solution...Appendix 2: The above code was saved on my hard disk under
~/stackexchange/answers_stackexchange_4625494.jland after making it executable withchmod +x(linux box, of course) it was started.Results:
Appendix 3:
I've implemented also the arithmetical approach, which comes in after a $(x,a,b)$ enumeration in some wanted range. The code was not successful, so it is skipped here.