ZSQRT := proc(a,u0,p) local u,pk,k,e,uk,i;
>     u := mods(u0,p);
>     i := modp(1/(2*u0),p);
>     e := a - u^2;
>     pk := p;
>     for k do
>         printf("e = %d\n", e);
>         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: