How fast can we find all Four-Square combinations that sum to N?

A question was asked at Stack Overflow here:

Given an integer $$N$$, print out all possible combinations of integer values of $$A,B,C$$ and $$D$$ which solve the equation $$A^2+B^2+C^2+D^2 = N$$.

This question is of course related to Bachet's Conjecture in number theory (sometimes called Lagrange's Four Square Theorem because of his proof). There are some papers that discuss how to find a single solution, but I have been unable to find anything that talks about how fast we can find all solutions for a particular $$N$$ (that is, all combinations, not all permutations).

I have been thinking about it quite a bit and it seems to me that it can be solved in $$O(N)$$ time and space, where $$N$$ is the desired sum. However, lacking any prior information on the subject, I am not sure if that is a significant claim on my part or just a trivial, obvious or already known result.

So, the question then is, how fast can we find all of the Four-Square Sums for a given $$N$$?

OK, here's the (nearly) $$O(N)$$ algorithm that I was thinking of. First two supporting functions, a nearest integer square root function:

// The nearest integer whose square is less than or equal to N
public int SquRt(int N) {
return (int) Math.Sqrt((double) N);
}

And a function to return all TwoSquare pairs summing from 0 to N:

// Returns a list of all sums of two squares less than or equal to N, in order.
public List<List<int[]>> TwoSquareSumsLessThan(int N) {
// Make the index array
List<int[]>[] Sum2Sqs = new List<int[]>[N + 1];

// Get the base square root, which is the maximum possible root value
int baseRt = SquRt(N);
for (int i = baseRt; i >= 0; i--) {
for (int j = 0; j <= i; j++) {
int sum = (i * i) + (j * j);
if (sum > N) {
break;
} else {
// Make the new pair
int[] sumPair = { i, j };
// Get the sumList entry
List<int[]> sumLst;
if (Sum2Sqs[sum] == null) {
// make it if we need to
sumLst = new List<int[]>();
Sum2Sqs[sum] = sumLst;
} else {
sumLst = Sum2Sqs[sum];
}

// Add the pair to the correct list
}
}
}

// Collapse the index array down to a sequential list
List<List<int[]>> result = new List<List<int[]>>();
for (int nn = 0; nn <= N; nn++) {
if (Sum2Sqs[nn] != null)
}

return result;
}

Finally, the algorithm itself:

// Return a list of all integer quads (a,b,c,d), where:
//      a^2 + b^2 + c^2 + d^2 = N,
// and  a >= b >= c >= d,
// and  a,b,c,d >= 0
public List<int[]> FindAllFourSquares(int N) {
// get all two-square sums <= N, in descending order
List<List<int[]>> Sqr2s = TwoSquareSumsLessThan(N);

// Cross the descending list of two-square sums <= N with
// the same list in ascending order, using a Merge-Match
// algorithm to find all combinations of pairs of two-square
// sums that add up to N
List<int[]> hiList, loList;
int[] hp, lp;
int hiSum, loSum;
List<int[]> results = new List<int[]>();
int prevHi = -1, prevLo = -1;

// Set the Merge sources to the highest and lowest entries in the list
int hi = Sqr2s.Count - 1;
int lo = 0;

//  Merge until done ..
while (hi >= lo) {
// check to see if the points have moved
if (hi != prevHi) {
hiList = Sqr2s[hi];
hp = hiList;     // these lists cannot be empty
hiSum = hp * hp + hp * hp;
prevHi = hi;
}

if (lo != prevLo) {
loList = Sqr2s[lo];
lp = loList;     // these lists cannot be empty
loSum = lp * lp + lp * lp;
prevLo = lo;
}

// Do the two entries' sums together add up to N?
if (hiSum + loSum == N) {
// they add up, so cross the two sum-lists over each other
foreach (int[] hiPair in hiList) {
foreach (int[] loPair in loList) {
// Make a new 4-tuple and fill it

// Only keep those cases where the tuple is already sorted
//(Otherwise it's a duplicate entry)
if (quad >= quad) { // (only need to check this one case, the others are implicit)
}

// (there's a special case where all values of the 4-tuple are equal
// that should be handled to prevent duplicate entries, but I'm
// skipping it for now)
}
}

// Both the HI and LO points must be moved after a Match
hi--;
lo++;
} else if (hiSum + loSum < N) {
lo++; // too low, so must increase the LO point
} else { // must be > N
hi--;   // too high, so must decrease the HI point
}
}

return results;
}

As I said before, it should be pretty close to $$O(N)$$. However, as Yuval Filmus points out, as the number of Four Square solutions to N can be of order $$N.ln(ln(N))$$, then this algorithim could not be less than that.

• Yes, please post it. I'm still developing the linear algorithm details, but I'm pretty sure it's valid. Aug 1 '12 at 21:04
• For the record, it appears that sometimes there are as many as $\Omega(N\log\log N)$ solutions, so we can't really have an $O(N)$ algorithm. Aug 2 '12 at 14:37
• From here, it looks like the catch (and the extra non-linear factor) comes from the two foreach() loops within your main while loop; your total time is basically $\displaystyle\sum_{i=0}^{N/2} |hiList_{N-i}|*|loList_i|$, and the problem is that the sizes of hiList and loList aren't necessarily bounded by any constant. Aug 2 '12 at 18:30
• Yes, that's correct, however your formula's a little bit off because first i ranges from 0 to apprx. NPI/8, and second only a fraction of the values of i satisfy hiList(N-i)+loList(i) = N, so they are not all added in. In any event, there's no way to fix this and I am pretty sure that this gives the minimum possible complexity of O(Nlog(log(N))). Aug 2 '12 at 21:11
• But we can have an algorithm that runs in O (max (N, "number of solutions")), taking O (n) space. Jan 13 '18 at 1:43

Juho's algorithm can be improved to an $O(N)$ algorithm using meet-in-the-middle. Go over all pairs $A,B \leq \sqrt{N}$; for each pair such that $M=A^2+B^2 \leq N$, store $(A,B)$ in some array $T$ of length $N$ (each position $M$ could contain several pairs, which might be stored in a linked list). Now go over pairs $M,N-M$ such that the corresponding cells in $T$ are non-empty.

This way we get an implicit representation of all quadruples. If we want to list all of them, then we can't do any better than $\Omega(N\log\log N)$, since Jacobi's four square theorem shows that (for odd $N$) the number of representations is $8\sigma(N)$, and there are infinitely many integers such that $\sigma(N) \geq (e^\gamma - \epsilon) N\log\log N$ (see Grönwall's theorem).

In order to get less trivial algorithms, one can try to factor $N$ over the appropriate quaternion ring, since we know that the representations as sums of two squares correspond (in some sense) to these factorizations, through Lagrange's four-square identity. We would still need to find all representations of any relevant prime.

• Hmm, the meet-in-the-middle thing sounds very similar to what I am working on (almost done) which is an ascending/descending Merge-Match algorithm over the TwoSquare pairs. Does that sound the same? Aug 2 '12 at 14:35
• It's probably the same, meet-in-the-middle is such a common heuristic that it must have many different names. Aug 2 '12 at 14:36
• Umm, I've been out of academia for thirty years, what's the $\sigma(N)$ thing mean? (or can you point me to a reference?) thnx. Aug 2 '12 at 14:40
• Or is that $\sigma(N)$ really a $\omicron(N)$? Aug 2 '12 at 15:01
• Sum of divisors function indeed. Aug 2 '12 at 15:48

The obvious algorithm that runs in quadratic time enumerates all tuples $$A,B,C,D \leq \sqrt[]{N}$$. This can be done in four loops, so the total time complexity becomes $$O(N^2)$$. It also clearly enumerates all solutions.

As relating algorithms, Rabin and Shallit  present two randomized algorithms for decomposing integers as sum of squares. For two squares, they give a $$O(\log^2 n)$$ expected time algorithm. For four squares, they give a $$O(\log^2 n \log \log n)$$ expected time algorithm. Note that the algorithms don't give you all the solutions, but merely just one.

• For a trivial algorithm, you only need loops for A, B, and C and then calculate D and check it is an integer. If you require A ≤ B ≤ C ≤ D you should get O (N^1.5) with a rather small constant. Jan 7 '17 at 13:12
• About 0.04 N^1.5 triples (A, B, C), and checking that N - A^2 - B^2 - C^2 is a square can be done very quickly. Jan 7 '17 at 22:36

The number of solutions is exactly $8\sum d$, where $d$ goes over all divisors of $n$ which are not multiples of 4. This is a theorem of Jacobi.