Four Squares Together Lagrange's four square theorem tells us any natural number can be represented as the sum of four square numbers. Your task is to write a program that does this.

Input: A natural number (below 1 billion)

Output: Four numbers whose squares sum to that number (order doesn't matter)

Note: You don't have to do a brute force search! Details here and here. If there  is a function that trivializes this problem (I will determine), it's not allowed. Automated prime functions and square root are allowed. If there is more than one representation any one is fine. If you chose to do brute force, it must run within reasonable time (3 minutes) 

sample input

123456789


sample output (either is fine)

10601 3328 2 0
10601 3328 2
Mathematica  61 66 51
Three methods are shown.  Only the first approach meets the time requirement.

1-FindInstance (51 char)
This returns a single solution the equation.
FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &amp;


Examples and timings
FindInstance[a^2 + b^2 + c^2 + d^2 == 123456789, {a, b, c, d}, Integers] // AbsoluteTiming


{0.003584, {{a -&gt; 2600, b -&gt; 378, c -&gt; 10468, d -&gt; 2641}}}

FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &amp;[805306368]



{0.004437, {{a -&gt; 16384, b -&gt; 16384, c -&gt; 16384, d -&gt; 0}}}


2-IntegerPartitions
This works also, but is too slow to meet the speed requirement.
f@n_ := Sqrt@IntegerPartitions[n, {4}, Range[0, Floor@Sqrt@n]^2, 1][[1]]

Range[0, Floor@Sqrt@n]^2 is the set of all squares less than the square root of n (the largest possible square in the partition).
{4} requires the integer partitions of n consist of 4 elements from the above-mentioned set of squares.
1, within the function IntegerPartitions returns the first solution.
[[1]] removes the outer braces; the solution was returned as a set of one element.

f[123456]


{348, 44, 20, 4}


3-PowerRepresentations
PowerRepresentations returns all of the solutions to the 4 squares problem.  It can also solve for sums of other powers.
PowersRepresentations returns, in under 5 seconds,  the 181 ways to express 123456789 as the sum of 4 squares:
n= 123456;
PowersRepresentations[n, 4, 2] //AbsoluteTiming


However, it is far too slow for other sums.