I've already received a lot of help with this problem (thank you all), and I've written up an answer with all of the information I've got in it as to why these approximations are so good.
My question is why the following rational function is such a good approximation for the square root of $x$?
$$\sqrt{x}\approx\frac{256+1792x+1120 x^{2}+112 x^{3}+x^{4}}{1024+1792x+448 x^{2}+16 x^{3}}$$
The following image illustrates that for $3\leqslant x\leqslant 7$ the two curves are almost identical:
More generally why does the following formula hold with such accuracy?:
$$\sqrt{n}\approx\frac{m^{8}+28 m^{6}n+70 m^{4}n^{2}+28 m^{2}n^{3}+n^{4}}{8 m^{7}+56 m^{5}n+56m^{3}n^{2}+8m n^{3}}$$
for any real number (not just integer) $n\geqslant 1$ and:
$$m=\left\lfloor\frac{1+\sqrt{4n-3}}{2}\right\rfloor$$
The following graph shows that for $0\leqslant n\leqslant10000$ the two curves cannot be separated by eye:
And this is the logarithm of the error in that approximation which decreases cyclically (presumably due to the floor in the definition of $m$) but then starts to become slightly less accurate stepwise as $n$ increases:
(Edit: Several people have mentioned Pade approximants below and that the first formula is a Pade approximant, which explains why it is such a good approximation. I would really especially like to know why the second formula is such a good approximation on the wide range I have shown in the second graph (is it related to Pade approximants also?), and I would really love to know why the process I outline in the section below of how I got these functions gives Pade approximants at all from ratios of binomial coefficients. At first glance, it seems as if there might be something deep going on here?)
How I found these and what I already know:
While messing around with maths once I noticed that $(3+\sqrt{8})^{8}\approx1331714$ with very good accuracy. Since the binomial theorem gives that $(3+\sqrt{8})^{8}=665857+235416\sqrt{8}$ and since $1331714=2\times665857$ I saw that that gave the approximation:
$$\sqrt{8}\approx\frac{665857}{235416}$$
which is quite accurate. I found this intriguing, but exploring further I found that similar near integer results were produced by:
$$(1+\sqrt{2})^{8}\approx1154$$
$$(2+\sqrt{3})^{8}\approx37634$$
$$(2+\sqrt{4})^{8}=65536$$
$$(2+\sqrt{5})^{8}\approx103682$$
$$(2+\sqrt{6})^{8}\approx153632$$
$$(3+\sqrt{7})^{8}\approx1032224$$
and so on, and I found that in general I could get a very accurate near integer from any (at least small) expression of the form $(m+\sqrt{n})^{8}$ where as $n$ increased over the integers $m$ increased as:
$$m=\left\lfloor\frac{1+\sqrt{4n-3}}{2}\right\rfloor$$
When these expressions were expanded out they were all in the form $a+b\sqrt{n}\approx2a$ where:
$$a=m^{8}+28 m^{6}n+70 m^{4}n^{2}+28 m^{2}n^{3}+n^{4}$$
$$b=8 m^{7}+56 m^{5}n+56m^{3}n^{2}+8m n^{3}$$
From this you get my original equations that I posted.
Now there is nothing particularly interesting about the eighth powers I have been working with; ninth powers work as well (slightly more accurately) and seventh powers are slightly less accurate, giving a similar formula:
$$\sqrt{n}\approx\frac{m^{7}+21 m^{5}n+35 m^{3} n^{2}+7m n^{3}}{7 m^{6}n+35 m^{4} n^{2}+21 m^{2} n^{3}+n^{4}}$$
Thus some of the reason behind these formulae must be the fact that the $m+\sqrt{n}$ have powers tending to integers, and just about the only thing I have been able to find out about why this might work is that some of these numbers are Pisot-Vijayaraghavan numbers (i.e. their powers tend towards integers).
However, this does not explain the phenomenon, since I do not know why the effect would be greatest for $m=\left\lfloor\frac{1+\sqrt{4n-3}}{2}\right\rfloor$ or importantly why the convergence should be of the form $(m+\sqrt{n})^{c}=a+b\sqrt{n}\approx 2a$ (which is where the approximation comes from) or why the formula derived should also hold for non-integral $n$. Rather the important part seems to lie in the formulation as:
$$\left(\frac{Even Binomial Terms}{Odd Binomial Terms}\right)^{\pm 1}$$
and perhaps in the expression $m=\left\lfloor\frac{1+\sqrt{4n-3}}{2}\right\rfloor$.
So my question is: why do these formulae work?
Best Answer
Concerning the first question $$\sqrt{x}\approx\frac{x^4+112 x^3+1120 x^2+1792 x+256}{16 x^3+448 x^2+1792 x+1024}$$ is the $\text{Pade}_{(4,3)}$ expansion of $\sqrt x$ centered at $x=4$ (have a look here).
Around a given point, any function can be approximated using Padé approximants which are ratios of polynomial of selected degrees.
If you consider $\text{Pade}_{(4,4)}$, it would be $$\sqrt{x}\approx \frac{9 x^4+336 x^3+2016 x^2+2304 x+256}{x^4+144 x^3+2016 x^2+5376 x+2304}$$ which is much better.
For the same number of coefficients, Padé approximants are incredibly more powerful than Taylor expansions. For example, in the case you gave, over the range $3\leq x \leq 5$, the largest error is less than $2.5 \times 10^{-9}$ that is to say more than ten times smaller than the result of a Taylor expansion to $O\left((x-4)^8\right)$.
We could play a lot with this kind of work and be sure that it is used a lot for many problems.
Edit (four years later)
Concerning the first part of the problem, we can make it more general and write (keeping the idea of a $[4,4]$ Padé approximant built at $x=a$) $$\sqrt[n] x=\sqrt[n] a\frac{1+b_1 t+b_2t^2+b_3t^3+b_4t^4}{1+c_1 t+c_2t^2+c_3t^3+c_4t^4} \quad \text{where} \quad t=\frac{x-a}{2 a n}$$ where the coefficients are $$b_1=4n+1$$ $$ b_2=\frac{3(3 n+1) (4 n+1)}{7}$$ $$b_3=\frac{2(2 n+1) (3 n+1) (4 n+1)}{21}$$ $$ b_4=\frac{(n+1) (2 n+1) (3 n+1) (4 n+1)}{105} $$ $$c_1=4n-1$$ $$ c_2=\frac{3(3 n-1) (4 n-1)}{7}$$ $$c_3=\frac{2(2 n-1) (3 n-1) (4 n-1)}{21}$$ $$ c_4=\frac{(n-1) (2 n-1) (3 n-1) (4 n-1)}{105} $$