How to find all possible values ​​of four variables with a quadratic sum up to N?

Question

A^2+B^2+C^2+D^2 = N Specify an integer N , print all possible combinations of integer values ABCD that solve the equation.

I guess we can do better than brute force.

Answer 1

There is interesting background information on the Wikipedia page, but the four-square Lagrange theorem (or rather, Bachet's theorem - only Lagrange proved it), does not go into details on how to find these squares.

As I said in my comment, the solution will not be trivial. This article discusses the decidability of four square sums. The document states that:

There is no convenient algorithm (with the exception of the simple one mentioned in the second paragraph of this document) for finding additional solutions that are indicated by calculating the representations, but perhaps this will simplify the search by giving an idea of ​​what types of solutions exist and do not exist.

There are some more interesting facts related to this topic. There are other theorems that claim that each integer can be written as the sum of four specific multiple squares. For example, each integer can be written as N = a ^ 2 + 2b ^ 2 + 4c ^ 2 + 14d ^ 2. There are 54 cases that are true for all integers, and Ramanujan provided a complete list in 1917.

See Modular Forms for more information. This is not easy to understand if you do not have any background in number theory. If you could generalize the forms of Ramanujan 54, it might be easier for you with this. With that said, in the first article I'm quoting, there is a small snippet that discusses an algorithm that can find any solution (although I find this to be a bit complicated):

For example, in 1911 it was reported that Gottfried Calculator Rell was asked to reduce N = 15663 as the sum of four squares. He produced a solution of 125 ^ 2 + 6 ^ 2 + 1 ^ 2 + 1 ^ 2 in 8 seconds, then immediately 125 ^ 2 + 5 ^ 2 + 3 ^ 2 + 2 ^ 2. A more complicated problem (reflected by the first term, which is farther from the original number, with correspondingly later terms) it took 56 seconds: 11399 = 105 ^ 2 + 15 ^ 2 + 8 ^ 2 + 5 ^ 2. In general, the strategy should begin with the first term being the largest square below N and try to present the smaller remainder as the sum of three squares. Then set the first term to the next largest square below N, etc. . Over time, a lightning calculator will become familiar with the expression of a small number as the sum of squares, which will speed up the process.

(Emphasize mine.)

The algorithm is described as recursive, but it can be easily implemented iteratively.

Answer 2

Naive brute force will look something like this:

 n = 3200724; lim = sqrt (n) + 1; for (a = 0; a <= lim; a++) for (b = 0; b <= lim; b++) for (c = 0; c <= lim; c++) for (d = 0; d <= lim; d++) if (a * a + b * b + c * c + d * d == n) printf ("%d %d %d %d\n", a, b, c, d); 

Unfortunately, this will lead to more than a trillion cycles, but not too efficient.

In fact, you can do it much better if you take into account the huge number of impossibilities at each level with something like:

 #include <stdio.h> int main(int argc, char *argv[]) { int n = atoi (argv[1]); int a, b, c, d, na, nb, nc, nd; int count = 0; for (a = 0, na = n; a * a <= na; a++) { for (b = 0, nb = na - a * a; b * b <= nb; b++) { for (c = 0, nc = nb - b * b; c * c <= nc; c++) { for (d = 0, nd = nc - c * c; d * d <= nd; d++) { if (d * d == nd) { printf ("%d %d %d %d\n", a, b, c, d); count++; } tot++; } } } } printf ("Found %d solutions\n", count); return 0; } 

This is still brute force, but not so cruel, because it understands when it is necessary to stop each level of the cycle as soon as possible.

In my (relatively) modest field, which takes the second ^(a) to get all the solutions for numbers up to 50,000. In addition, it starts to take more time:

  n time taken ---------- ---------- 100,000 3.7s 1,000,000 6m, 18.7s 

For n = ten million , it lasted about an hour and a half before I killed him.

So, I would say that brute force is quite acceptable to a certain point. In addition, more mathematical solutions will be required.

---

For greater efficiency, you can only check those solutions where d >= c >= b >= a . This is because you could collect all the solutions from these combinations into permutations (with possible duplication, where the values ​​of two or more of a , b , c or d identical).

In addition, the body of cycle d does not need to check each value of d , only the last possible.

Getting results for 1,000,000 in this case takes less than ten seconds, and no more than six minutes:

  0 0 0 1000 0 0 280 960 0 0 352 936 0 0 600 800 0 24 640 768 : : : : 424 512 512 544 428 460 500 596 432 440 480 624 436 476 532 548 444 468 468 604 448 464 520 560 452 452 476 604 452 484 484 572 500 500 500 500 Found 1302 solutions real 0m9.517s user 0m9.505s sys 0m0.012s 

This code follows:

 #include <stdio.h> int main(int argc, char *argv[]) { int n = atoi (argv[1]); int a, b, c, d, na, nb, nc, nd; int count = 0; for (a = 0, na = n; a * a <= na; a++) { for (b = a, nb = na - a * a; b * b <= nb; b++) { for (c = b, nc = nb - b * b; c * c <= nc; c++) { for (d = c, nd = nc - c * c; d * d < nd; d++); if (d * d == nd) { printf ("%4d %4d %4d %4d\n", a, b, c, d); count++; } } } } printf ("Found %d solutions\n", count); return 0; } 

---

And, according to DSM's suggestion, cycle d can disappear altogether (since there is only one possible value of d (discounting negative numbers), and it can be calculated), which leads to one million cases up to two seconds for me, and ten million cases for a much more manageable 68 seconds.

This version is as follows:

 #include <stdio.h> #include <math.h> int main(int argc, char *argv[]) { int n = atoi (argv[1]); int a, b, c, d, na, nb, nc, nd; int count = 0; for (a = 0, na = n; a * a <= na; a++) { for (b = a, nb = na - a * a; b * b <= nb; b++) { for (c = b, nc = nb - b * b; c * c <= nc; c++) { nd = nc - c * c; d = sqrt (nd); if (d * d == nd) { printf ("%d %d %d %d\n", a, b, c, d); count++; } } } } printf ("Found %d solutions\n", count); return 0; } 

---

^(a) : All timings are done with internal printf commented out, so I / O does not distort the numbers.

Answer 3

It seems that all integers can be made using this combination:

 0 = 0^2 + 0^2 + 0^2 + 0^2 1 = 1^2 + 0^2 + 0^2 + 0^2 2 = 1^2 + 1^2 + 0^2 + 0^2 3 = 1^2 + 1^2 + 1^2 + 0^2 4 = 2^2 + 0^2 + 0^2 + 0^2, 1^2 + 1^2 + 1^2 + 1^2 + 1^2 5 = 2^2 + 1^2 + 0^2 + 0^2 6 = 2^2 + 1^2 + 1^2 + 0^2 7 = 2^2 + 1^2 + 1^2 + 1^2 8 = 2^2 + 2^2 + 0^2 + 0^2 9 = 3^2 + 0^2 + 0^2 + 0^2, 2^2 + 2^2 + 1^2 + 0^2 10 = 3^2 + 1^2 + 0^2 + 0^2, 2^2 + 2^2 + 1^2 + 1^2 11 = 3^2 + 1^2 + 1^2 + 0^2 12 = 3^2 + 1^2 + 1^2 + 1^2, 2^2 + 2^2 + 2^2 + 0^2 . . . 

etc.

As I started working in my head, I thought that these would be only perfect squares that have more than one possible solution. However, after listing them, it seems to me that there is no obvious order. However, I thought about which algorithm I consider most suitable for this situation:

It is important to use a 4-tuple (a, b, c, d). In any given 4-tuple, which is the solution a ^ 2 + b ^ 2 + c ^ 2 + d ^ 2 = n, we establish the restriction that a is always the largest of 4, b is the following, and so on and so forth ., as:

 a >= b >= c >= d 

We also note that a ^ 2 cannot be less than n / 4, otherwise the sum of the squares must be less than n.

Then the algorithm:

 1a. Obtain floor(square_root(n)) # this is the maximum value of a - call it max_a 1b. Obtain the first value of a such that a^2 >= n/4 - call it min_a 2. For a in a range (min_a, max_a) 

At this moment, we chose a special a and now look at bridging the gap from a ^ 2 to n - ie (n - a ^ 2)

 3. Repeat steps 1a through 2 to select a value of b. This time instead of finding floor(square_root(n)) we find floor(square_root(n - a^2)) 

etc. etc. Thus, the whole algorithm will look something like this:

 1a. Obtain floor(square_root(n)) # this is the maximum value of a - call it max_a 1b. Obtain the first value of a such that a^2 >= n/4 - call it min_a 2. For a in a range (min_a, max_a) 3a. Obtain floor(square_root(n - a^2)) 3b. Obtain the first value of b such that b^2 >= (n - a^2)/3 4. For b in a range (min_b, max_b) 5a. Obtain floor(square_root(n - a^2 - b^2)) 5b. Obtain the first value of b such that b^2 >= (n - a^2 - b^2)/2 6. For c in a range (min_c, max_c) 7. We now look at (n - a^2 - b^2 - c^2). If its square root is an integer, this is d. Otherwise, this tuple will not form a solution 

In steps 3b and 5b, I use (n - a ^ 2) / 3, (n - a ^ 2 - b ^ 2) / 2. We divide by 3 or 2, respectively, because the number of values ​​in the tuple is not yet “fixed”.

Example:

doing this with n = 12:

 1a. max_a = 3 1b. min_a = 2 2. for a in range(2, 3): use a = 2 3a. we now look at (12 - 2^2) = 8 max_b = 2 3b. min_b = 2 4. b must be 2 5a. we now look at (12 - 2^2 - 2^2) = 4 max_c = 2 5b. min_c = 2 6. c must be 2 7. (n - a^2 - b^2 - c^2) = 0, hence d = 0 so a possible tuple is (2, 2, 2, 0) 2. use a = 3 3a. we now look at (12 - 3^2) = 3 max_b = 1 3b. min_b = 1 4. b must be 1 5a. we now look at (12 - 3^2 - 1^2) = 2 max_c = 1 5b. min_c = 1 6. c must be 1 7. (n - a^2 - b^2 - c^2) = 1, hence d = 1 so a possible tuple is (3, 1, 1, 1) 

These are the only two possible sets: hey presto!

Answer 4

Nebff has a great answer. one suggestion:

  step 3a: max_b = min(a, floor(square_root(n - a^2))) // since b <= a 

max_c and max_d can also be improved.

Here is another attempt:

 1. generate array S: {0, 1, 2^2, 3^2,.... nr^2} where nr = floor(square_root(N)). 

Now the task is to find 4 numbers from the array that sum (a, b, c, d) = N;

 2. according to neffa post (step 1a & 1b), a (which is the largest among all 4 numbers) is between [nr/2 .. nr]. 

We can loop a from nr to nr / 2 and compute r = N - S [a]; Now the question is to find 3 numbers from S: the sum (b, c, d) = r = N -S [a];

here is the code:

 nr = square_root(N); S = {0, 1, 2^2, 3^2, 4^2,.... nr^2}; for (a = nr; a >= nr/2; a--) { r = N - S[a]; // it is now a 3SUM problem for(b = a; b >= 0; b--) { r1 = r - S[b]; if (r1 < 0) continue; if (r1 > N/2) // because (a^2 + b^2) >= (c^2 + d^2) break; for (c = 0, d = b; c <= d;) { sum = S[c] + S[d]; if (sum == r1) { print a, b, c, d; c++; d--; } else if (sum < r1) c++; else d--; } } } 

runtime O(sqare_root(N)^3) .

Here is the test result that runs java on my virtual machine (time in milliseconds, result # - total number of valid combination, time 1 with print, time2 without printout):

 N result# time1 time2 ----------- -------- -------- ----------- 1,000,000 1302 859 281 10,000,000 6262 16109 7938 100,000,000 30912 442469 344359