Why does this square root algorithm work? I've been doing some elliptic curve cryptography, and a library I'm using has this slightly bizarre algorithm for computing modular square roots:

Let $x$ be some quadratic residue modulo $p$ for some large prime $p$. Let $h$ be some element of $\mathbb{Z} \setminus p\mathbb{Z}$ such that $h^2-4x$ has no square root modulo $p$. Let $V_n$ be the Lucas sequence with $V_0=2$, $V_1=h$ and $V_n=hV_{n-1}+xV_{n-2}$. Let $k = \frac{p+1}{2}$. Then $V_k = 2\sqrt{x}$. This is obviously less efficient than the well known algorithm for $p \cong 3 \mod 4$, but apparently works in a more general case

Immediately, one might believe this definition to be contradictory; indeed the characteristic equation of this recurrence is solved by square rooting $h^2-4x$, but, being a Lucas sequence, this impossible term never appears in any particular value of $V_k$.

The only reference to this algorithm I've found is as an exercise in an undergraduate-level textbook, which suggests that this is solvable by some fairly simple maths, but nothing I've tried seems to result in a solution -- computing $V_{\frac{p+1}{2}}^2$ yields $V_{\frac{p+1}{2}}^2 = V_2 + 2x^{\frac{p+1}{2}} = h^2 - 2x + 2x^{\frac{p+1}{2}}$, which looks almost nothing like the $4x$ I would expect.
This is Cipolla's algorithm, I believe.

The integers mod $p$ we call $\mathbf{F}_p$, and since $h^2 -4x$ is a quadratic nonresidue mod $p$, the polynomial $P(Y) := Y^2-(h^2-4x) \in \mathbf{F}_p[Y]$ has no roots. We can then mod out the polynomial ring $\mathbf{F}_p[Y]$ to get $\mathbf{K} := \mathbf{F}_p[Y]/(Y^2-(h^2-4x))$ (which by some algebra we know is isomorphic to $\mathbf{F}_{p^2}$). It is in this field $\mathbf{K}$ that $h^2-4x$ has a square root (one can think of it as the indeterminate $Y = \sqrt{h^2-4x}$)

In this extension field $\mathbf{K}$ (which is still characteristic $p$, so $(m+n)^p = m^p+n^p$ for all $m,n\in \mathbf{K}$) we have that $\left(h + \sqrt{h^2-4x}\right)^p = h^p + (\sqrt{h^2-4x})^p$. 

By Euler's criterion, $(h^2-4x)^{\frac{p-1}{2}} \cong -1 \pmod p$ and Fermat's little theorem $h^p \cong h \pmod p$, and all field operations in $\mathbf{F}_p$ are the same in $\mathbf{K}$ on $\mathbf{F}_p$ elements. Therefore  $h^p + (\sqrt{h^2-4x})^p = h + (h^2-4x)^{\frac{p-1}{2}} \sqrt{h^2-4x} = h - \sqrt{h^2-4x}$. 

By now, with all that math we conclude that $\left(h + \sqrt{h^2-4x}\right)^p = h - \sqrt{h^2-4x}$. The same argument also yields $\left(h - \sqrt{h^2-4x}\right)^p = h + \sqrt{h^2-4x}$. 

Finally we get to the square root of $x$. Clearly $$4x = h^2 - (h^2 - 4x) = (h + \sqrt{h^2-4x})(h - \sqrt{h^2-4x}) = (h+\sqrt{h^2-4x})(h+\sqrt{h^2-4x})^p = (h+\sqrt{h^2-4x})^{p+1}$$

Now, the $k$th term of the Lucas sequence $V_k := hV_{k-1} + xV_{k-2}$ is $V_k = a^k + b^k$ where $a = \frac{h + \sqrt{h^2-4x}}{2}$ and $b = \frac{h - \sqrt{h^2-4x}}{2}$. So $$V_{\frac{p+1}{2}} = \left(\frac{h + \sqrt{h^2-4x}}{2}\right)^{\frac{p+1}{2}} + \left(\frac{h - \sqrt{h^2-4x}}{2}\right)^{\frac{p+1}{2}}$$

$$V_{\frac{p+1}{2}} = \left(\frac{1}{2}\right)^{\frac{p+1}{2}} \left( (h + \sqrt{h^2-4x})^{\frac{p+1}{2}} + (h - \sqrt{h^2-4x})^{\frac{p+1}{2}} \right)$$

$$V_{\frac{p+1}{2}} = \left(\frac{1}{2}\right)^{\frac{p+1}{2}} ( \sqrt{4x} + \sqrt{4x}) = \left(\frac{2}{2^{\frac{p-1}{2}}}\right) \sqrt{x} = \pm 2 \sqrt{x}$$