1st Pari/GP experience, determine sum of squares for prime p = 1 (mod 4), implement "assert()"

sq2.gp

\\ sq2(p) determines sum of squares for prime p = 1 (mod 4), implement "assert()"
\\
\\ based on
\\ - https://pari.math.u-bordeaux.fr/Events/PARI2019b/talks/init.pdf
\\ - https://pari.math.u-bordeaux.fr/Events/PARI2019b/talks/prog.pdf
\\ - https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.15.1/tutorial.pdf
\\ - https://skalatan.de/en/archive/pariguide/doc/Programming_under_GP.html
\\ - https://home.gwu.edu/~maxal/gpscripts/binomod.gp
\\
assert(b, v, s) = { if(!(b), error(Str(v) Str(s))); }

bits(N) = { #binary(N); }
digits_(N) = { #digits(N); }

\\ from Robert Chapman Python code (https://math.stackexchange.com/a/5883/1084297)
\\
mods(a,n) =
{
   assert(n>0, n, " <=0");
   a = a % n;
   if(2*a>n, return(a-n));
   return(a);
}

powmods(a,r,n) =
{
   my(out = 1);
   while(r>0,
     if(r%2==1, 
       r-=1;
       out = mods(out*a, n);
     );
     r = floor(r/2);
     a = mods(a*a, n);
   );
   out;
}

quos(a, n) =
{
   assert(n>0, n, " <= 0");
   floor((a - mods(a, n))/n);
}

grem(w,z) =
{
   my(w0,w1,z0,z1,n);
   [w0,w1] = w;
   [z0,z1] = z;
   n = z0*z0 + z1*z1;
   assert(n>0, "division by ", n);
   u0 = quos(w0*z0+w1*z1,n);
   u1 = quos(w1*z0-w0*z1,n);
   [w0-z0*u0+z1*u1, w1-z0*u1-z1*u0];
}

ggcd(w,z) =
{
   while(z!=[0,0],
     [w,z] = [z, grem(w,z)];
   );
   w;
}

root4m1(p) =
{
   my(k,j,a,b);
   assert(p>1 && p%4==1, p, " not 1 (mod 4)");
   k=floor(p/4);
   j=2;
   while(1,
     a=powmods(j,k,p);
     b=mods(a*a,p);
     if(b==-1, return(a););
     assert(b==1, p, " not prime");
     j+=1;
   );
}

sq2(p) =
{
   my(a,x,y);
   assert(p>1 && p%4==1, p, " not 1 (mod 4)");
   a=root4m1(p);
   [x,y] = ggcd([p,0], [a,1]);
   [abs(x), abs(y)];
}

sq2(p) determines sum of squares for prime p = 1 (mod 4).
Calling with 97 works (97 == 9^2 + 4^2), with non-prime 85 asserts (with assert() from line 8):

sq2() is quite fast, here on prime factors of RSA-768 ("##" reports time taken by previous command):

? sq2(33478071698956898786044169848212690817704794983713768568912431388982883793878002287614711652531743087737814467999489)
[5777499308315937806714612090100083182548223876592695640833, 313964076552969675277689223019594269926166116787107977840]
? ##
  ***   last result: cpu time 4 ms, real time 6 ms.
?
? sq2(36746043666799590428244633799627952632279158164343087642676032283815739666511279233373417143396810270092798736308917)
[4045487632634994038632464260708513408124495115537369857049, 4514429474584457604059260685088937606791715198364890782254]
? ##
  ***   last result: cpu time 4 ms, real time 4 ms.
? 

Revision 2 corrects function definitions, adds important GP programming links.