Phil Carmody <firstname.lastname@example.org> wrote: > Pubkeybreaker <email@example.com> writes: >> On Aug 11, 7:27?am, David Bernier <david...@videotron.ca> wrote: >>> >>> Concerning Fermat's sum of two squares theorem, that every >>> prime p of the form p = 4n+1 can be expressed as >>> the sum of two perfect squares, books and on-line >>> references don't say much about the computational >>> complexity of decomposing a generic prime of the form >>> 4n+1 as a sum of two squares. >>> >>> The PARI-gp package finds that the 51-digit number N51: >>> N51 = 314,159,265,358,979,323,846,264,338,327,950,288,419,716,939,937,673 >>> >>> is a prime p of the form 4n+1 >>> >>> Is it feasible today to have a decent chance >>> of decomposing a prime such as N51 as a sum of two squares ? >> >> Yes. It is trivial. Look up "Cornacchia's Algorithm" > > Sometimes Cornacchia-Smith. It seems as if Smith may have > had the idea earlier, but they are as near as dammit. > > In pari, the cheat or cheap method to decompose > such a prime P into 2 squares is to factor P+0*I. > The answer will come back (a+b*I)*(a-b*I), so a^2+b^2=P. > I have no idea what algorithm it uses, so it may not > be efficient asymptotically, but I've never seen it > break a sweat (OK, it does a probable primality test > first, but that's not so expensive).
This can be done simply by using the Euclidean algoritm in Z[i] to compute gcd(p,i-j) = a+bi, where j^2 = -1 (mod p). Then p = a^2+b^2. Below is a Macsyma program that does such. See also my posts here http://math.stackexchange.com/questions/5877
/* For p=1 (mod 4) compute w = a+bi with p = ww' = a^2 + b^2 */
psplit(p) := block([j:0, e:floor(p/4), algebraic:true, x, y ],
block([modulus:p], /* search for j = sqrt(-1) (mod p) */ for a:2 while j^2 # -1 do j:rat(a)^e ),
x:p, y : j - %i, /* find gcd(p,j-i) in Z[i] by Euclidean algorithm */
do ( if y = 0 then return(x), j : rat(x - y * floor(rat(x/y))), x:y, y:j ))$