# Efficient computation of integer representation as 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.

-
For a fast probabilistic approach, based on Cramer's conjecture, simply take small values of $z$ in the appropriate arithmetic progression, and check if $n-z^2$ is a prime number which is $1$ mod $4$. Then, to get $p=x^2+y^2$, use Lagrange's algorithm. – Dror Speiser Aug 9 '12 at 6:59
Sounds reasonable, but I have a couple of questions. First, notice that your proposition is true only for $m \equiv 1 \pmod{4}$. In that case $x$ and $z$ are odd, $y$ is even, and therefore $m - z^2 = x^2 + y^2 \equiv 1 \pmod{4}$. Is it always possible to find such a prime $p$ so that $p = x^2 + y^2$? Do you know any proof for this? Second, for $m \equiv 3 \pmod{8}$ we have $m - z^2 = x^2 + y^2 \equiv 2 \pmod{8}$. So in this case, I guess, we should look for $p$ such that $2p = x^2 + y^2$. Again, why does this $p$ always exists? Also, could you please give a link to an article about algorithm? – Anton Aug 9 '12 at 7:33
I believe there are no references in the literature to this problem, nor does the current knowledge on primes allow to prove theorems such as those above. I also have a copy of Grosswald's book, and I can assure you that it does not contain an answer to this problem. I hate to sound so pessimistic, but hey, on the bright side, if you do find anything better than exponential complexity, you'll probably be the first to do so! – Dror Speiser Aug 9 '12 at 8:51
@Dror factoring is subexponential, so finding a $z$ that does work will give subexponential complexity or am I missing something? This approach is far from efficient though. – joro Aug 9 '12 at 9:06
You are missing two things: the first is that what you said isn't exactly true, you don't know how long you'll search for a $z$ that actually works, i.e. that $n-z^2$ is a sum of two squares - sure, once it works you have $x,y$ from factorisation, but how many $z$ must you try before one works? second, subexponential factoring are probabilistic algorithms, without proofs, just like the method that has been proposed. – Dror Speiser Aug 9 '12 at 9:15

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.

-
Looks like knowing one solution allows full? parametrization over Q. – joro Oct 22 '12 at 12:44

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]);
);
}

-
OK, this works. The question is, does this work in a polynomial time? There exist two algorithms (by Serret and by Schoof) that compute $x$, $y$ for a given $p=x^2+y^2$ in a polynomial time. But how many values of z will we check before finding a proper one? I need this for cryptographic purposes, so efficiency is crucial. – Anton Aug 9 '12 at 8:52
By the way, here's the online tool for solving the two squares problem: numbertheory.org/php/main_pell.html. Also, I found a book that contains an algorithm by Serret that solves the two squares problem in $O(\log^4{p})$ (based on ERH): plouffe.fr/simon/math/…. – Anton Aug 9 '12 at 8:57
The algorithm is probabilistic and I don't know how many values of z are needed. You can check the program for your data (probably lowering 10^8 to 10^6, check for yourself). As the program showed 1332 bit integer was solved in about a minute. – joro Aug 9 '12 at 9:00
That sounds like a pretty good result. Ok, I'll do some more research on how fast this works and how many $z$ should be checked and write about it later. By the way, I think that it is more efficient to start not with $z=0$, but with $z= \lfloor n / \sqrt{3} \rfloor$. Then we'll get the smallest value of $n−z^2$ and it would be more probable for this number to be prime. – Anton Aug 9 '12 at 9:14
Here's some heuristics for predicting how many values of $z$ are needed: the solutions to $x^2 + y^2 + z^2 = n$ are roughly equidistributed on the sphere, and the number of solutions is proportional to the class number $h(-4n)$, which grows like $\sqrt{n}$. So consider $C \cdot \sqrt{n}$ points equidistributed on the sphere of area $4 \pi n$. What is the statistical expectation of the smallest (positive) $z$-coordinate among these points? – Marty Aug 9 '12 at 22:07