I found a very neat ancient hindi formula for approximating square roots using rational numbers. After doing some algebra on the formula, i came across with this recursive relation:
Given any number $k\in\ \mathbb{Q}$: $$\sqrt{k}≈x_{n+1}=\frac{x_n^4+6kx_n^2+k^2}{4(x_n^3+kx_n)}$$ In the original article, it is stated that successive iterations of this relation result in a cuartic convergence, so this formula is very sharp. However, for this formula to actually work, we need to make an initial guess: We have to pick some $x_0\in \mathbb{Q}$ such that for a sufficiently small open interval…let’s say $I=(\frac{\sqrt{|k|}}{2};\sqrt{|k|})\subset \mathbb{Q}$: $$x_0\in I$$ This is great, because now we can have a way to generalize root approximation for some algebraic numbers. But the question arises: is there any way to actually “compute” $x_0$ algebraically? The best description i came off with is the floor function $x_0=\lfloor\sqrt{|k|}\rfloor$ But that would not be an actual algebraic expression for this, as i would like. Can you come up with anything better?
Looking at the linked paper and at your formula, it seems to me that this is almost the same as the original Householder formula.
Let me write it in a different manner : we want to find the zero of function $$f(x)=x^2-k$$ and we shall approximate it as its $[1,m]$ Padé approximant $P_m$ around a given value $x=a$. Notice that this is the same as a Newton-like method of order $\color{red}{(m+2)}$.
This write $$f(x) \sim P_m=\frac {(a^2-k) + \alpha_m (x-a)}{1+\sum_{p=1}^m b_p\,(x-a)^p }$$ which gives $$x_{(m)}=a+\frac {k-a^2}{\alpha_m}$$
I shall write the first cases making the result as a single fraction as you did (notice that $a$ corresponds to your $x_n$ and $x_{(m)}$ to your $x_n$)
$$x_{(0)}=\frac{a^2+k}{2 a} \tag 0$$ $$x_{(1)}=\frac{a(a^2+3k) }{3 a^2+k}\tag 1$$ $$x_{(2)}=\frac{a^4+6 a^2 k+k^2}{4 a \left(a^2+k\right)}\tag 2$$ which is your formula.
But we can continue and generate $$x_{(3)}=\frac{a \left(a^4+10 a^2 k+5 k^2\right)}{5 a^4+10 a^2 k+k^2}\tag 3$$ $$x_{(4)}=\frac{a^6+15 a^4 k+15 a^2 k^2+k^3}{a(6 a^4+20 a^2 k+6 k^2)}\tag 4$$ $$x_{(5)}=\frac{a \left(a^6+21 a^4 k+35 a^2 k^2+7 k^3\right)}{7 a^6+35 a^4 k+21 a^2 k^2+k^3}\tag 5$$ $$x_{(6)}=\frac{a^8+28 a^6 k+70 a^4 k^2+28 a^2 k^3+k^4}{8 a \left(a^6+7 a^4 k+7 a^2 k^2+k^3\right)}\tag 6$$ All of the above are just particular cases of Householder general approach which, symbolically, write $$x_{n+1} = x_n + d\; \frac {\left(1/f\right)^{(d-1)}(x_n)} {\left(1/f\right)^{(d)}(x_n)}$$
Notice that is $k$ and $a$ are rational numbers, all iterates are also rational.
Quoting the linked paper
"In the examples presented in the Bakhshali manuscript, this algorithm is used to obtain rational approximations to square roots only for integer arguments $q$, only for integer-valued starting values $x_0$, and is only applied once in each case"
With regard to the starting value, using Darboux theorem for Newton method, since $f''(x) >0$, we shall avoid overshoot of the solution if the first guess $a$ is such that $f(a) >0$ that is to say $a^2 > k$.
Notice that I asked here what about the higher order methods; this question (which is four years old) never received any comment or answer.