Explanation/Prove for this Continued fraction Algorithm Computing regular continued fractions by iteratively inverting the remainder is easiest, but runs into problems with fixed precision arithmetics. Its possible to work with arbitrary precision libraries, but that is really slow.
Due to other posts on this page I came along the algorithm described here. It is very intuitive, but it requires a repeated computation of square roots (or rather the integer part of the square root), which still seems unnecessary:
I also found the algorithm hidden on this page, which reads implemented in Python:
def cf_sqrt(D):
    a0 = int(sqrt(D))
    result = [a0]

    an, Pn, Qn = a0, 0, 1
    while an != 2*a0:
        Pn = an*Qn - Pn
        Qn = (D - Pn**2)/Qn
        an = int((a0 + Pn)/Qn)
        result.append(an)
    return result

It only needs a single square root evaluation and other than that only basic arithmetic operations. However, I can not figure out, why this actually works. I can verify the result for individual numbers, but I'd like to have a prove, that this functions indeed computes the continued fraction of $\sqrt{D}$.
 A: This is not really an answer, but let me illustrate a fairly standard procedure for getting the expansion of $\sqrt n$ where $n>0$, nonsquare. I’ll do it for $\sqrt7$.
Preliminary step, find $m=\lfloor\sqrt7\rfloor$, i.e. the largest integer for which $m^2<7$, so $m=2$.
Then compute:
\begin{align}
\sqrt7-2=\frac{\sqrt7-2}1&=\frac3{\sqrt7+2}&=\frac1{(\sqrt7+2)/3}&=\frac1{1+(\sqrt7-1)/3}\\
\frac{\sqrt7-1}3&=\frac2{\sqrt7+1} &=\frac1{(\sqrt7+1)/2}&=\frac1{1+(\sqrt7-1)/2}\\
\frac{\sqrt7-1}2&=\frac3{\sqrt7+1}&=\frac1{(\sqrt7+1)/3}&=\frac1{1+(\sqrt7-2)/3}\\
\frac{\sqrt7-2}3&=\frac1{\sqrt7+2}&=\frac1{4+\sqrt7-2}\\
\sqrt7-2=\text{( it repeats )}
\end{align}
Thus the expansion of $\sqrt7$ is $2+\frac1{1+}\frac1{1+}\frac1{1+}\frac1{4+\cdots}$ . When I said, “precise arithmetic, derationalizing the denominator when necessary”, I meant exactly the step that takes, for instance, $(\sqrt7-1)/3$ to $(7-1)/\bigl(3(\sqrt7+1)\bigr)=2/(\sqrt7+1)$. It’s just the process you learned in high-school to “rationalize the denominator”, except going backwards.
I’m pretty sure that if you analyse the process step by step, you’ll see that it’s performed exactly by the Python routine that you’ve quoted.
