From: rusin@vesuvius.math.niu.edu (Dave Rusin)
Newsgroups: sci.math
Subject: Re: A problem
Date: 22 Oct 1997 18:51:40 GMT

In article <344B39FB.A212A244@iol.it>, piSuKE  <pisuke@iol.it> wrote:
>x,y are positive integers such that x^2+y^2=1997. Find x,y!

"By inspection", {x,y}={29,34}

>How do I solve x^2+y^2=p , where p is prime and 4 divides (p-1)?

You can easily do this in polynomial time. (That is, the number of steps
is a polynomial in  log(p).)

1. Compute a square root  r  of (-1) mod p. If you select an  x  at
random from {1, 2, ..., p-1}, it will be a nonsquare mod p with 
probability 1/2; you will know when you compute  x^((p-1)/2), since this
must be  +1  for a square and -1 for a nonsquare. Your chances of
getting a square in  log_2(p)  choices of  x  are only  1/p, which is
to say, nearly impossible for large  p. So it's relatively easy to
find a nonsquare  x. Then  r=x^((p-1)^2)  is a square root of (-1).
(Note that computing  x^m  mod  p  takes at most  O(log(m)) multiplies mod p).

2. Now you know  r^2+1 = N * p  for some  N. Now, the  x  and  y  you
seek are essentially unique, since they are the integers for which
x + i y  is a prime divisor of  p  in the ring of Gaussian integers  Z[i].
Thus the equation to which you have come  may be written
	(r+i)*(r-i) = N*(x+iy)*(x-iy)
Using the fact that  x+iy  is prime (and the fact that  Z[i]  is a Euclidean
domain!) we deduce that  x+iy  divides one of  r+i  and  r-i; equivalently,
either  x+iy  or x-iy  divides  r+i.  Therefore, one of these two primes
is a common divisor of  r+i  and  p. It must then divide the GCD of there
two Gaussian integers. (In fact,   p  being a product of just two Gaussian
primes, not both dividing  r+i,  the GCD  must be _equal to_  x+iy.)

3. So now use the Euclidean algorithm to compute the GCD of  p  and  r+i:
gcd(a+bi, c+di):=gcd(c+di, (a+bi)-(c+di)*(x+iy)), where  x+iy is the
nearest Gaussian integer to  
	(a+bi)/(c+di) = (ac-bd)/(c^2+d^2) + i (ad+bc)/(a^2+b^2)


Here's Maple code:

> p:=157526249:
> 2 &^ ((p-1)/2) mod p;  #yields  1  --  2's a square.
> 3 &^ ((p-1)/2) mod p;  #yields  -1  --  3's not a square.
> r:=3&^((p-1)/4) mod p; #gives the sqrt of -1.
                                 r := 76631264
> ggcd:=proc(a,b) local x:
> if b=0 then a: else x:=round(Re(a/b))+I*round(Im(a/b)):ggcd(b,a-b*x):fi:
> end:
>
> ggcd(p,r+I);
                                 9532 + 8165 I

> 9532^2+8165^2;
                                   157526249

I just learned this over the summer (when I really Had To Know); thanks to
Jeff Shallit for leading me to this (obvious!) technique.

dave