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

Question

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
                sumLst.Add(sumPair);
            }
        }
    }

    // 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)
            result.Add(Sum2Sqs[nn]);
    }

    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[0];     // these lists cannot be empty
            hiSum = hp[0] * hp[0] + hp[1] * hp[1];
            prevHi = hi;
        }

        if (lo != prevLo) {
            loList = Sqr2s[lo];
            lp = loList[0];     // these lists cannot be empty
            loSum = lp[0] * lp[0] + lp[1] * lp[1];
            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
                    int[] quad = new int[4];
                    quad[0] = hiPair[0];
                    quad[1] = hiPair[1];
                    quad[2] = loPair[0];
                    quad[3] = loPair[1];

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

                    // (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.

Answer 1

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.

Answer 2

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 [1] 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.

---

[1] M. O. Rabin, J. O. Shallit, Randomized Algorithms in Number Theory, Communications on Pure and Applied Mathematics 39 (1986), no. S1, pp. S239–S256.

Answer 3

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.