RCS: I'm not sure about the rigor, but a very practical method of finding 3-square reps is to look at N-k^2, with k decrementing from floor(sqrtN)). Remove small prime factors, give up when there are an odd number of 3s or 7s, and hope that the cofactor is a 4K+1 prime. WDS: actually, I was just going to suggest exactly that. Now it seems plausible that the chance a number Q is prime is of order 1/lnQ hence this would find about sqrt(N)/ln(N) representations, and with the mean time between reps-found being of order logN. However, that reasoning was not rigorous. But it dawns on me that we can in fact get both rigor and decent speed as follows: ALGORITHM: Choose k at random in {1,2,...,floor(sqrt(N))}. If N-k^2 happens to be representable as a sum of 2 squares, you are done. Otherwise try again. It is a theorem of C.L.Siegel, re-proved in 1 page by Dorian M. Goldfeld: A simple proof of Siegel's theorem, Proc. Nat. Acad. Sci. U.S.A. 71 (1974) 1055 http://www.pnas.org/content/71/4/1055.full.pdf+html that the "class number" is lower bounded by constant*N^(1/2-epsilon) which tells us as a corollary (via work by Gauss) that every number N which is representable as a sum of 3 squares has at least constant*N^(1/2-epsilon) representations for any fixed epsilon>0. We also know you can factor N into primes in N^o(1) randomized time. We also know that numbers N have at most N^o(1) representations as sum of 2 squares. Putting all those things together, we conclude THEOREM: Above algorithm will run in expected N^o(1) time to find a representation of N as sum of 3 squares (if it has any). (Also, N is representable as the sum of 3 squares iff N is NOT of the form 4^a * (8*k+7) so it is trivial to decide representability.) Also Erich Hecke proved [proof published by E.Landau: Gottingen Nachr. (1918) 285-295] assuming the extended Riemann hypothesis that the number of 3-square reps was >constant*sqrt(N)/logN. So we might hope algorithm above will run in O(logN) expected primality tests, except that perhaps for a finite set of N it won't work. That's not a theorem, it's merely a hope, at present. RCS: If you have one rep as A^2+B^2+C^2, you can walk to more: a 5step move is see if any of the pairs A^2+B^2, A^2+C^2, or B^2+C^2 is divisible by 5; if so, multiply by (3+-4i)/5 to find the other 2-square rep for that pair; rinse & repeat. A 3step move is more interesting: fudge the signs so that 3 | (A+B+C); then subtract 2*(A+B+C)/3 from each of A,B,C. There are 7steps, etc.; roughly a commutative group; usually one or two kinds of step will cover all the 3-square representations. 3-step example: 1001 = 30^2 + 10^2 + 1^2; signs +30+10-1; 39*(2/3)=26; (30,10,-1) - (26,26,26) = (4,-16,-27); 16+256+729 = 1001. Gauss showed relationship between 3-square reps and binary quadratic forms; the steps are equivalent to composition of BQFs by constant forms. Steps exist for more general ternary quadratic forms like A^2+A*B+2*B^2+3*C^2. WDS: excellent. RCS: also said stuff about continued fractions and gaussian integers vis-a-vis the 4-square problem... WDS: indeed, for that see G. Rousseau: On a construction for the representation of a positive integer as the sum of four squares, L'Enseignement Math. (2) 33,3-4 (1987) 301-306.