This algorithm due to Lagrange gives the expansion as a continued fraction of a quadratic number
(u+√D)/v, with u, v, D integers, D not a square.

Working

The algorithm requires v dividing D-u².
If v doesn't divide D-u², we can replace u by u|v|, D by Dv² and v by v|v| then v|v| divides Dv²-u²v²,
although being the same quadratic number.

Let U_{0}=u, V_{0}=v, a=E[√D],
With E[x]=integer part of x (ceil function).
And added variables X, Y, P X_{-1}=-u, X_{0}=v
Y_{-1}=1, Y_{0}=0
P_{-1}=0, P_{0}=1

until sequence V_{n} becomes periodic.
The a_{i} are the terms of the continued fraction,
P_{i}/Y_{i} are the convergents
and X_{i}²-DY_{i}²=(-1)^{i}V_{0}V_{i}