Question: Hoxw to decompose a prime number into the sum of 2 squares

I am trying to decompose an isprime into sum of 2 squares.
Can you tell me why yhse procedure are not goog.

Sumof2Squares:= proc(p::And(prime, satisfies(p-> irem(p,4)=1)))
local x, y:= 1;
   x:= mods(Roots(x^2+y^2), p)[2,1];
   while x^2+y^2 > p do
      (x,y):= FermatDescent(x,y,p)
   end do;
end proc:

FermatDescent:= proc(x::posint, y::posint, p::posint)
   m:= (x^2+y^2)/p,
   a:= mods(x,m),  
   b:= mods(y,m);

   (abs((a*x+b*y)/m), abs((a*y-b*x)/m))
end proc:

Thank you.

Please Wait...