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

Problem Detail: 

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;         int 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.

Asked By : RBarryYoung

Answered By : Yuval Filmus

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.

Best Answer from StackOverflow

Question Source : http://cs.stackexchange.com/questions/2988