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}$.