ZSQRT := proc(a,u0,p) local u,pk,k,e,uk,i;
>     u := mods(u0,p); uk := u;
>     i := modp(1/(2*u0),p);
>     e := a - u^2;
>     pk := p;
>     for k do
>         printf("got %d = %d expected\n", e, a-u^2);
>         if e = 0 then return(u); fi;
>         if e < 0 then return(FAIL) fi;
>         uk := mods( iquo(e,pk)*i, p );
>         u := u + uk*pk;
> 
>         e := e + (-2*u*uk*pk-uk^2*pk^2);
>         
>         pk := p*pk;
>     od;
> end: