Computer Science Stack Exchange is a question and answer site for students, researchers and practitioners of computer science. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
Yes, please post it. I'm still developing the linear algorithm details, but I'm pretty sure it's valid. – RBarryYoung 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. – Yuval Filmus 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. – Steven Stadnicki 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))). – RBarryYoung Aug 2 '12 at 21:11
up vote 13 down vote accepted

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.

share|cite|improve this answer
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? – RBarryYoung 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. – Yuval Filmus 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. – RBarryYoung Aug 2 '12 at 14:40
Or is that $\sigma(N)$ really a $\omicron(N)$? – RBarryYoung Aug 2 '12 at 15:01
Sum of divisors function indeed. – Yuval Filmus Aug 2 '12 at 15:48

I think a $o(N^2)$ time algorithm is not a trivial one and requires some insight if one exists. 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.

share|cite|improve this answer

The number of solutions is exactely 8* sum d, where d is a divisor of n and 4 does not divide d. This is a theorem of Jacobi.'s_four-square_theorem

share|cite|improve this answer
And how does this answer the question? The task is to give all these quadruples! – Raphael Aug 19 '15 at 20:47
This is already mentioned in my answer. – Yuval Filmus Aug 19 '15 at 21:40

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.