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) {
            } 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[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.

  • $\begingroup$ Yes, please post it. I'm still developing the linear algorithm details, but I'm pretty sure it's valid. $\endgroup$ Aug 1, 2012 at 21:04
  • 5
    $\begingroup$ 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. $\endgroup$ Aug 2, 2012 at 14:37
  • 1
    $\begingroup$ 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. $\endgroup$ Aug 2, 2012 at 18:30
  • $\begingroup$ 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))). $\endgroup$ Aug 2, 2012 at 21:11
  • $\begingroup$ But we can have an algorithm that runs in O (max (N, "number of solutions")), taking O (n) space. $\endgroup$
    – gnasher729
    Jan 13, 2018 at 1:43

3 Answers 3


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.

  • $\begingroup$ 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? $\endgroup$ Aug 2, 2012 at 14:35
  • 1
    $\begingroup$ It's probably the same, meet-in-the-middle is such a common heuristic that it must have many different names. $\endgroup$ Aug 2, 2012 at 14:36
  • $\begingroup$ 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. $\endgroup$ Aug 2, 2012 at 14:40
  • $\begingroup$ Or is that $\sigma(N)$ really a $\omicron(N)$? $\endgroup$ Aug 2, 2012 at 15:01
  • 1
    $\begingroup$ Sum of divisors function indeed. $\endgroup$ Aug 2, 2012 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 [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.

  • $\begingroup$ 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. $\endgroup$
    – gnasher729
    Jan 7, 2017 at 13:12
  • $\begingroup$ 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. $\endgroup$
    – gnasher729
    Jan 7, 2017 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.

  • 1
    $\begingroup$ And how does this answer the question? The task is to give all these quadruples! $\endgroup$
    – Raphael
    Aug 19, 2015 at 20:47
  • 1
    $\begingroup$ This is already mentioned in my answer. $\endgroup$ Aug 19, 2015 at 21:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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