Efficient computation of integer representation as a sum of three squares Recently I've been studying the problem of integer representation as sum of three squares. Most of the articles that I've found study the function $r_m(n)$ which counts the number of representations of $n$ as the sum of $m$ squares. However, this is not what I am interested in. What I'm looking for is an efficient way (for some given $n$) to find $x$, $y$ and $z$ such that $n = x^2 + y^2 + z^2$. I need to find at least one such representation. Can you recommend me some articles that study this problem?
P.S. I believe that Emil Grosswald's book "Representation of integers as Sums of Squares" contains the answer. However, I could not find this book on my university's web-site.
 A: This problem is discussed in my paper with Rabin, Randomized algorithms in number theory,
Commun. Pure Appl. Math. 39, 1985, S239 - S256.  We give an algorithm that, assuming a couple of reasonable conjectures, will produce a representation as a sum of three squares in polynomial time.
A: A modification of Dror's comment.
This probabilistic algorithm worked for me.
The main idea is to pick some $z$, compute $m=n - z^2$, factor $m$ with trial division and express it as a sum of two squares if possible. The probability of finding prime $m=4a+1$ or $2p$ is high enough for practical purposes.
The algorithm:


*

*z:=0

*z:=z+1

*m:=n-z^2

*if can't trial factor m goto 2

*if m=x^2+y^2 (the factorization is known) then x^2+y^2+z^2=n. Done

*goto 2


Added Experimental results:  100 random integers of size 1000 bits were solved at average time 3.5 seconds per solution. 10 random integers of size 10000 bits were solved at average time 5.8 minutes.
Here is a pari/gp program and example:
/*
? n=nextprime(10^220+30000)*nextprime(2^1000+40000000);n%8
%136 = 3
? t=threesquares(n);

? ##
  ***   last result computed in 2,210 ms.
? round(log(n))
%138 = 1200
*/

pl=10^6;\\ bound for trial division, may need change
default(primelimit,pl);
{
twosquares(n)=
local(K,i,v,p,c1,c2);
K=bnfinit(x^2+1);
v=bnfisintnorm(K,n);
for(i=1,#v,p=v[i];c1=polcoeff(p,0);c2=polcoeff(p,1);if(denominator(c1)==1&&denominator(c2)==1,return([c1,c2])) );
return([]);
}

{
threesquares(n)=
local(m,z,i,x1,y1,j,fa,g);
if(n/4^valuation(n,4)==7,return([]););
for(z=1,n,
\\forstep(z=sqrtint(n),1,-1,
m=n-z^2;
if(m%4==3,next);
print1(z," ",);
fa=factor(m,pl);
g=1;
for(i=1,#fa~,if(!ispseudoprime(fa[i,1])/*||!isprime(fa[i,1])*/||(fa[i,2]%2==1&&fa[i,1]%4==3),g=0;break; ));
if(!g,next);
print("\nfound ",z," "," m=",m,factor(m));
j=twosquares(m);
print("j=",j);
x1=abs(j[1]);
y1=abs(j[2]);
return([x1,y1,z]);
);
}

