Re: how to best determine a^2+b^2=-1 (mod p) for odd prime p?

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Mon, Apr 27, 2026 at 04:52:46PM +0200, hermann@stamm-wilbrandt.de wrote:
> I learned here long ago how to determine a^2+b^2=1 (mod p) for prime
> p=4*k+1:
> 
> sosp1(p)={
>   [M,V]=halfgcd(lift(sqrt(Mod(-1,p))),p);
>   return([V[2],M[2, 1]]);
> }
> 
> Today I heard in a seminar about existence of a^2+b^2=-1 (mod p) for any odd
> prime p.

> Asking gemini it said what I implemented, but I used sqrt(Mod(-1-x^2,p)) for
> Mod(-1-x^2,p) that is a quadratic residue instead of "Tonelli Shanks"
> (perhaps that is used in GP sqrt?).

Yes, sqrt uses a variant of Tonelli-Shanks.

> sosm1(p)={
>   while(1,x=random(p);w=Mod(-1-x^2,p);if(kronecker(lift(w),p)==1,
>     return([x,lift(sqrt(w))])));
> }
> 
> Works, and since half of numbers generated by random(p) are quadratic
> residues of p,
> expected number of loop executions is 2.
> 
> Question is, whether there is a GP builtin function to compute a,b?

No, because there are too many solutions!

An alternative solution (which is probably slower) for p=3 mod 4 but
does not require sqrt:

sosm1(p)=
{
  my(r,a=ffgen((t^2+1)*Mod(1,p)));
  until(norm(r)==-1,
    r=random(a)^((p-1)/2));
  Vec(r.pol);
}

or even simply:

sosm1(p)=
{
  my(r);
  until(norm(r)==-1,
    r=((random(p)+I*random(p))*Mod(1,p))^((p-1)/2));
  lift([real(r),imag(r)]);
}

Cheers,
Bill