Theorem of Lagrange: If N>0 is an integer then it has a representation as a sum of 4 squares N=a^2+b^2+c^2+d^2. Theorem of Euler: If, further, N's squarefree part is divisible only by primes 2 and 4*k+1, then (and only then) N=a^2+b^2. But: How can one QUICKLY FIND a,b,c,d given N?. Answer: Using GCDs over the Gaussian integers a+i*b and the quaternions a+i*b+j*c+k*d respectively! Here are the details. N=SUM OF TWO SQUARES: If N is square we are done. If N=2 then 2=1^2+1^2. If N is a prime of form 4*k+1 then use Shanks-Tonelli algorithm to solve x^2+1=0 (mod N) for x. Then find the Gaussian integer GCD(x+i, N). The result is a+i*b where a^2+b^2=N. Otherwise factor N into primes and squares, apply above methods to represent all the primes and squares, and multiply all the Gaussian integers you get. N=SUM OF FOUR SQUARES: We may factor N to find its squarefree part. So wlog we need only consider squarefree N. But actually I think it will suffice merely to remove all factors of 4 from N (which is way faster than factoring N to remove *all* squares) in the below. And now: If N=odd, then 2N-1+4*k*N is (for some k) a prime that is 1 mod 4, and hence of form a^2+b^2 where a and b found as above. If N=even, then, since N=squarefree, N=2(mod 4), in which case N-1+4*k*N is prime and 1 mod 4 (for some k), hence is expressible as a^2+b^2 as above. [Actually it is not necessary to choose k to get a prime since some composite N that are 1 mod 4 also are expressible as sum of two squares. But primes are frequent so that random k's will yield primes often, and primality testing is fast.] Let g=a+b*i. We know by construction that N divides a^2+b^2+1. Now take the quaternion (right) GCD(N, g+j) to find a representation of N as the sum of 4 squares. These algorithms run in time polynomial in the number of bits in N provided: * the prime factorization of N is available (for the 2 squares problem) * we are allowed to use a random number generator (for the 4 squares problem). What about 3 squares?