Kunerth's algorithm from 1878, for determining modular sqrt

Kunerth.gp

assert(b,s)=if(!(b), error(Str(s)));

m=eval(getenv("m"));
b=m.mod;
c=lift(m);

if(ispseudoprime(c),s=sqrt(Mod(-b,c)),issquare(Mod(b,-c),&s));
V=r=lift(s);;
print("V=",r);

e=simplify(((c*z+r)^2+b)/c);
f=1;
w=polcoeff(e,0);
if(type(w)!="t_INT"||w<0,e=simplify(((c*z+r)^2-b)/c);f=-1;w=polcoeff(e,0));
print("e=",e);
print("f=",f);
assert(issquare(w,&W));
print("W=",W);

beta=-(V\W);
alpha=W*(V+W*beta);
print("alpha=",alpha);
print("beta=",beta);

xx=alpha^2*x^2+(2*alpha*beta-f*b)*x+(beta^2-c);
print("xx=",xx);
nfr=nfroots(,xx);
print("nfr=",nfr);

{
  foreach(nfr,X,
    Y=Mod(alpha*X+beta,b);
    if(lift(Y^2)==c,
      print("X=",X);
      print("Y=",Y);
      print("Y^2=",Y^2)));
}

print("all asserts OK");

write("/dev/stderr", "\nb=m.mod; c=lift(m); V=r=lift(\"sqrt\"(Mod(-b,c)));");
write("/dev/stderr", "e=simplify(((c*z+r)^2±b)/c); f=±1; W=sqrt(polcoeff(e,0));");
write("/dev/stderr", "beta=-(V\\W); alpha=W*(V+W*beta);");
write("/dev/stderr", "xx=alpha^2*x^2+(2*alpha*beta-f*b)*x+(beta^2-c);  nfr=nfroots(,xx);");
write("/dev/stderr", "\n∀X∈nfr: Y=Mod(alpha*X+beta,b);");

Thanks for recovery this nice algorithm which is quite forgotten these days.

I've tried

$ m="Mod(861,1189)" gp -q < Kunerth.gp 
V=0
e=861*z^2 - 29/21
f=-1
  ***   at top-level: assert(issquare(w,&W))
  ***                 ^----------------------
  ***   in function assert: if(!(b),error(Str(s)))
  ***                               ^--------------
  ***   user error: 0
W=W
alpha=0
beta=0
xx=1189*x - 861
  ***   at top-level: nfr=nfroots(,xx)
  ***                     ^------------
  *** nfroots: incorrect type in nfrootsQ [not in Z[X]] (t_POL).
nfr=nfr
  ***   at top-level: foreach(nfr,X,Y=Mod(alpha*X+beta,b);if(lift(Y^
  ***                 ^----------------------------------------------
  *** foreach: incorrect type in foreach (t_POL).
all asserts OK

b=m.mod; c=lift(m); V=r=lift("sqrt"(Mod(-b,c)));
e=simplify(((c*z+r)^2±b)/c); f=±1; W=sqrt(polcoeff(e,0));
beta=-(V\W); alpha=W*(V+W*beta);
xx=alpha^2*x^2+(2*alpha*beta-f*b)*x+(beta^2-c);  nfr=nfroots(,xx);

∀X∈nfr: Y=Mod(alpha*X+beta,b);

and 123^2 = 861 (mod 1189)