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

[hermann@stamm-wilbrandt.de]

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?).

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?
I did not find some in GP doc.

Regards,

Hermann.