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:

sqrt-around-half

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):

sqrt-iterations

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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: