Efficient computation of integer representation as a sum of three squares

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 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:

1. z:=0
2. z:=z+1
3. m:=n-z^2
4. if can't trial factor m goto 2
5. if m=x^2+y^2 (the factorization is known) then x^2+y^2+z^2=n. Done
6. 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]);
);
}