Suppose $$$p$$$ is an odd prime and $$$a$$$ is a quadratic residue modulo $$$p$$$. Cipolla's algorithm shows that $$$b:=x^{(p+1)/2}\bmod{(x^2-tx+a)}$$$ such that $$$b^2=a$$$ for some irreducible polynomial $$$x^2-tx+a\in\mathbb{F}_p[x]$$$.
I don't know how to prove this properly. If there are any mistakes or typos, please let me know, thanks! Let $$$\alpha$$$ be a zero of $$$x^2-tx+a$$$, $$$\mathbb{F}_p(\alpha):=\lbrace a_0+a_1\alpha :a_0,a_1\in\mathbb{F}_p\rbrace$$$, and let $$$\beta$$$ be a zero of $$$x^2-(t^2-4a)$$$, $$$\mathbb{F}_p(\beta):=\lbrace a_0+a_1\beta :a_0,a_1\in\mathbb{F}_p\rbrace$$$. We may easily find homomorphisms between $$$\mathbb{F}_p(\alpha)$$$ and $$$\mathbb{F}_p(\beta)$$$. So I think we just need to show that $$$((t+\beta )/2)^{p+1}=a$$$.
A lot of blogs show that $$$(t+\beta)^p=t-\beta$$$ according to binomial theorem and Fermat's little theorem, so
Bostan and Mori's paper shows that the computation of $$$x^n$$$ modulo a monic polynomial sometimes is equivalent to the computation of one term of a rational function $$$P(x)/Q(x)$$$ where $$$P,Q$$$ are both polynomials and $$$\deg(P(x))\lt \deg(Q(x))$$$.
We have
and
This may save some multiplications.
Actually, the initial value of $$$k_1$$$ does no matter to the result. So we could mainly focus on the computation of one term of the inverse of formal power series. I don't know if we could do better.
Auto comment: topic has been updated by hly1204 (previous revision, new revision, compare).
Berlekamp-Rabin algorithm looks simpler, easier to implement and faster to me. Also it's easily generalized to finding roots of higher order and even roots of arbitrary polynomials over $$$\mathbb Z_p$$$.
Problem. For prime $$$p$$$, find $$$a$$$ such that $$$a^2 \equiv y \pmod p$$$.
Solution. Calculate $$$(z+x)^{\frac{p-1}{2}} \pmod{x^2-y}$$$ for random $$$z \in \mathbb Z_p$$$ such that $$$z^2 \neq y$$$.
If the result is $$$a_0 + a_1 x$$$ where $$$a_1 \neq 0$$$, then $$$a_0=0$$$ and $$$a=\frac{1}{a_1}$$$ is the solution.
Tl;dr.
Proof by AC.
Explanation.
Due to Chinese remainder theorem, it is equivalent to computing $$$(z+x)^{\frac{p-1}{2}}$$$ modulo $$$x-a$$$ and $$$x+a$$$.
Modulo $$$x-a$$$ it would be equal to $$$(z+a)^{\frac{p-1}{2}}$$$ and modulo $$$x+a$$$ it would be equal to $$$(z-a)^{\frac{p-1}{2}}$$$.
If both are $$$1$$$ or both are $$$-1$$$, the result would be $$$1$$$ or $$$-1$$$ modulo $$$x^2-y$$$ as well.
However if $$$z+a$$$ is a quadratic residue while $$$z-a$$$ is not (hence one makes $$$1$$$, the other makes $$$-1$$$), then
Actual probability of success is at least $$$\frac{1}{2}$$$, you can find a proof in some book or paper or whatever.
btw, subject can be shortened to this and be pretty fast, but ofc it is impossible to reproduce from scratch