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);");

Kunerth's algorithm is described here:
https://en.wikipedia.org/wiki/Kunerth%27s_algorithm

1$ ebook "Adolf Kunerth’s Modular Square Root Algorithm Explained: with examples in Mathematica":
https://www.amazon.com/Kunerths-Modular-Square-Algorithm-Explained-ebook/dp/B09K3TCQR1/ref=sr_1_1

German language pages 327-346 from 1878 in this comment:
https://math.stackexchange.com/questions/95443/how-to-compute-modular-square-roots-when-modulus-is-non-prime#comment9832085_252186
Adolf Kunerth, "Sitzungsberichte. Academie Der Wissenschaften" vol 78(2), 1878,

• p 327-338 (for quadratic equation algorithm),
• pp. 338–346 (for modular quadratic algorithm), available at Ernest Mayr Library, Harvard University

Kunerth.gp was created by porting Mathematica code fragments from the ebook to PARI/GP:
https://pari.math.u-bordeaux.fr/archives/pari-users-2311/msg00057.html

Example 2 from ebook (examples 3 and 4 from ebook in above posting).
The modulus 3872 is non-prime and does not get factored in calculating sqrt(353) (mod 3872).
Kunerth.gp is not a general modular square root determination algorithm.
Explanations after "all asserts OK" are sent to stderr; if you don't want them, append " 2>/dev/null".

$ m="Mod(353,3872)" gp -q < Kunerth.gp
V=94
e=353*z^2 + 188*z + 36
f=1
W=6
alpha=24
beta=-15
xx=[x - 8, 1; 36*x + 1, 1]
X=8
Y=Mod(177, 3872)
Y^2=Mod(353, 3872)
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=factor(alpha^2*x^2+(2*alpha*beta-f*b)*x+(beta^2-c));

i∈{1,2}:  X=-polcoeff(xx[i,1],0])/polcoeff(xx[i,1],1); Y=Mod(alpha*X+beta,b);
$  

Here is a new example calculation.
f=-1 is forced for getting a t_INT constant term of e, by subtracting b in simplify() instead of adding.
Modified f value influences computation of xx:

$ m="Mod(60,5*17)" gp -q < Kunerth.gp
V=25
e=60*z^2 + 50*z + 9
f=-1
W=3
alpha=3
beta=-8
xx=[x + 4, 1; 9*x + 1, 1]
X=-4
Y=Mod(65, 85)
Y^2=Mod(60, 85)
X=-1/9
Y=Mod(20, 85)
Y^2=Mod(60, 85)
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=factor(alpha^2*x^2+(2*alpha*beta-f*b)*x+(beta^2-c));

i∈{1,2}:  X=-polcoeff(xx[i,1],0])/polcoeff(xx[i,1],1); Y=Mod(alpha*X+beta,b);
$ 

Constant term of e was square 9, 5*17-60=25 was square as well.

Here is another square difference resulting in square constant of e and results:

$ m="Mod(13*17-7^2,13*17)" gp -q < Kunerth.gp 2>/dev/null
V=93
e=172*z^2 + 186*z + 49
f=-1
W=7
alpha=14
beta=-13
xx=[4*x - 3, 1; 49*x + 1, 1]
X=3/4
Y=Mod(108, 221)
Y^2=Mod(172, 221)
X=-1/49
Y=Mod(113, 221)
Y^2=Mod(172, 221)
all asserts OK
$ 

Bill Allombert made me aware that I can use nfroots(,P) to find the rational roots of P.
That simplified my code to determine the roots "by hand" significantly.
Diff for revisions 5 and 6 shows the changes, "nfroots()" made "factor()" call for xx unnecessary, and simplified code with foreach().
New variable nfr containing results of nfroots() call is printed as well:

$ m="Mod(60,85)" gp -q < Kunerth.gp
V=25
e=60*z^2 + 50*z + 9
f=-1
W=3
alpha=3
beta=-8
xx=9*x^2 + 37*x + 4
nfr=[-4, -1/9]~
X=-4
Y=Mod(65, 85)
Y^2=Mod(60, 85)
X=-1/9
Y=Mod(20, 85)
Y^2=Mod(60, 85)
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);
$