Exact randomized modular square roots in expected polynomial time

Reduce $a$ modulo $p$. If $p=2$, return its residue; if $a=0$, return zero. Otherwise [F1] rejects exactly the nonsquares by testing $a^{(p-1)\mathrm{/}2}=-1$. In the remaining case write $p-1=2^{s}q$ with $q$ odd and $s\ge 1$. Repeatedly sample $z$ uniformly from $1,\dots ,p-1$ until $z^{(p-1)\mathrm{/}2}=-1$. Uniform sampling is implemented by drawing $\lceil {\log }_{2}p\rceil$ bits and rejecting zero and integers at least $p$. Each draw is accepted with probability at least $1\mathrm{/}4$, and each accepted $z$ is a nonsquare with probability $1\mathrm{/}2$ by [F2]; thus this search has bounded expected trials and terminates almost surely.

Set $x=a^{(q+1)\mathrm{/}2}$, $t=a^q$, $c=z^q$ and $M=s$, all modulo $p$. Then $x^2=at$, $c$ has exact multiplicative order $2^M$, and $t^{2^{M-1}}=1$ by the tests in step 1.1. While $t\mathrm{\ne }1$, find the least $m\ge 1$ with $t^{2^m}=1$ by successive squaring; it satisfies $m<M$. Put $b=c^{2^{M-m-1}}$ and replace $(x,t,c,M)$ by $(xb,t{b}^2,b^2,m)$. The identity $x^2=at$ is preserved. The old $t$ has order exactly $2^m$, so $t^{2^{m-1}}=-1$: in a field, the only roots of $X^2-1$ are $1,-1$ by factorization. Likewise $b^2$ has order $2^m$ and its $2^{m-1}$-st power is $-1$. Therefore the new $t$ has order dividing $2^{m-1}$, while the new $c$ has order $2^m$. The invariants are maintained and $M$ strictly decreases. After at most $s$ iterations, $t=1$ and $x^2=a$. Return $x$.

Put $B=1+\lceil {\log }_{2}p\rceil$. Every modular integer has $O(B)$ bits. Schoolbook multiplication and division use $O(B^2)$ bit operations; binary repeated squaring computes a power with $O(B)$-bit exponent in $O(B)$ modular multiplications. These bounds follow directly by forming at most $B$ shifted rows for multiplication and eliminating at most $B$ quotient digits for division. Step 1.1 uses an expected bounded number of such powers. Each of the at most $B$ iterations in step 2.1 uses at most $O(B)$ modular multiplications to find $m$ and compute $b$. Its total cost is $O(B^4)$. Initial reduction has cost polynomial in the original length of $a$ and $p$. Thus the algorithm has the asserted expected bit bound; only the nonsquare search is randomized, and its outcome never changes the correctness of the final answer.