## Yet Another Square Root Algorithm (part III)

Finding good, fast, and stable numerical algorithms is usually hard, even for simple problems like the square root, subject we’ve visited a number of times. One method we haven’t explored yet, is using the Taylor series. The Taylor series for $\sqrt{n}$ around $a$ is given by:

$\displaystyle \sqrt{n}=\sqrt{a}+\frac{1}{2}\frac{n-a}{\sqrt{a}}-\frac{1}{8}\frac{(n-a)^2}{a^{3/2}}+\frac{1}{16}\frac{(n-a)^3}{a^{5/2}}+\cdots$

Can we exploit that to compute efficiently, and with some precision, square roots?

Well, yes. Kind of. I’ll come back later on the rather important caveat. But let’s try to find some regularity in the expansion. The coefficients seems cumbersome: $1, \frac{1}{2}, -\frac{1}{8}, \frac{1}{16}, -\frac{5}{128}, \frac{7}{256}, -\frac{21}{1024}, \frac{33}{2048}, \ldots$ until we remember that we’ve already saw them before. This allows us to write:

$\displaystyle \sqrt{n}=\sqrt{a}\sum_{i=0}^\infty \binom{\frac{1}{2}}{i}\left(\frac{n-a}{a}\right)^i$,

which is both regular and (relatively) simple. We can even factor the binomial coefficients and get an efficient iterative algorithm. Here’s a simple Mathematica implementation:

R[a_, n_, k_] := Module[{i, s = 0, c = Sqrt[a]},
For[i = 0, i <= k, i++,
s += c;
c *= 1/( i + 1) (1/2 - i) ((n - a)/a );
];
Return[s]
]


*
* *

Let’s see how this thing behaves. Let’s start with numbers on the interval $]0,1]$, with $a=\frac{1}{2}$. With three terms, we get:

with the green line being the true square root and the red dashed line being the approximation using the Taylor series. Not bad, but we see that the approximation gets somewhat rapidly worse as we get away from one half. That’s not unexpected: we used the expansion around $a=\frac{1}{2}$, and we use only three terms. When we are clearly far from $a=\frac{1}{2}$, the approximation is just plain wrong, and we must use a lot more terms and another $a$.

When we move away, we must chose different $a$ for different $n$. The following graph shows different $a$ giving different approximations (with different number of terms):

The graph gives a false sense of precision. The solutions are good to only a (very) few decimals.

*
* *

Choosing $a$ is a problem. The convergence radius is small, and a single $a$ can’t be used for a wide range of $n$. Worst, the Taylor expansion seems to ask us to compute $\sqrt{a}$, which is kind of what we’re trying to do in the first place. One way out of this problem is to use a $a$ that is a known perfect square, of about the same magnitude than $n$.

*
* *

Maybe we could use the first few terms to get a rough estimate for $\sqrt{n}$ and use that estimate to kick-start some other algorithm, say, Newton’s, to get a refined solution.