3

I need to compute the square roots of lots of numbers. The numbers increase monotonically by fixed step. For example, 1, 2, 3, ..., 1 000 000.

What is the fastest way to do so? Is it possible somehow to take the advantage of the growth and calculate the square roots incrementally?

I thought I would calculate the derivative and add it at each step. But the derivative contains again a square root so it doesn't help. Also, I don't mind to trade some precision for speed.

Ecir Hana
  • 1,025
  • 2
    But the derivative contains exactly the square root that you just calculated? In fact all derivatives can be calculated from that square root with just basic arithmetic, so you can use a higher-order expansion if you want. There's also the sequence $x_{n+1}=(x_n+a/x_n)/2$, which converges very quickly; you can use the previous square root as a good initial value; the advantage is that errors don't accumulate as in the other approach. But standard square root functions are already very well optimized; it will be hard to beat them with any approach like this. – joriki Jan 26 '13 at 10:12
  • @joriki: I tried this but I cannot get it to work, I always get slightly different result than the one from proper sqrt implementation. Please, do you know why? Here is what I do: http://codepad.org/UcwDA8ux – Ecir Hana Jan 26 '13 at 11:23
  • I wouldn't have expected that to work. You're using just a first-order approximation to estimate the square roots, and letting the errors accumulate without any corrections. I'm surprised it works as well as it does. I wasn't saying that something like this should work; I was just disagreeing with your statement that the derivative doesn't help because you'd need the square root for it. But you won't get enough precision just from the first derivative; you need to either use higher derivatives, or, as in my answer, correct through iteration. – joriki Jan 26 '13 at 11:29
  • @joriki ah, ok! Thanks! – Ecir Hana Jan 26 '13 at 11:31
  • @joriki btw, is there any difference when saying "higher derivates" and "forward differencing"? Or do they mean the same thing? – Ecir Hana Jan 26 '13 at 18:15
  • +1 for that sequence. Since it's based off of Newton's method, it has quadratic convergence. If you want better accuracy, just use more terms in the $x_{n+1}=(x_n+a/x_n)/2$ series. I suspect that setting the initial guess to the previously computed value will make this at least competitive with standard square root functions, especially as your sequence gets to very large numbers. – KevinG Jan 26 '13 at 19:53

4 Answers4

5

I'd suggest to calculate the square root $s_k\stackrel!=\sqrt{a_k}$ by using $s_{k-1}+(s_{k-1}-s_{k-2})$ as an initial value for the iteration $x_{n+1}=(x_n+a_k/x_n)/2$.

joriki
  • 243,700
  • 2
    Sorry to bother you, but could you explain to me what the symbol $\stackrel{!}{=}$ means please? – Thomas Russell Jan 26 '13 at 10:41
  • 2
    @Shaktal: It means (or at least I'm using it in that sense) "should be equal to", i.e. it's not an equality that's given but one that you're trying to achieve in what follows. – joriki Jan 26 '13 at 11:23
2

Here's a neat way to compute the next square root once the numbers are large. Let $s_n=\sqrt{n}$ be the previous value you computed and then use the Taylor expansion $$ s_{n+1}=\sqrt{n}\sqrt{1+1/n}=s_n\sqrt{1+1/n}=s_n\left(1+\frac{1}{2n}-\frac{1}{8n^2}+\frac{1}{16n^3}-\frac{5}{128n^4}+\cdots\right) $$ For example for $n=100$ the first three terms give you six correct digits. For binary arithmetic note the convenient powers of two in the denominators. The error will accumulate from one term to the next, but if we correct the answer in the cases $s_{m^2}=m$, then the error can't get too out of hand. Also note that anytime $n$ is not a prime, say $n=kl$ you can of course use $s_n=s_ks_l$ to compute the square root. I would at least do this for all even numbers, and then perhaps the odd numbered ones could be computed using the formula above.

I have another idea for how to make it faster using a rational approximation instead of a Taylor approximation, but I will need some time to work out the details.

Peder
  • 2,252
  • An improvement to the above suggestion is to compute $s_{n+1}$ from $s_{n-1}$ instead. We have $$s_{n+1}s_{n-1}=n\sqrt{1-1/n^2}=n\left(1-\frac{1}{2n^2}-\frac{1}{8n^4}-\frac{1}{16n^6}-\cdots\right)$$ As before the even numbers are computed from the odd ones using $\sqrt{2n}=\sqrt{2}\sqrt{n}$ iteratively until $n$ is odd. – Peder Jan 27 '13 at 17:18
1

Another approach: approximate the square root function by a polynomial, or maybe a rational function. The degree will depend on what accuracy you need and over how wide an interval. Calculating a good approximation is not easy, but you only have to do it once. If any of your numbers are close to zero, a polynomial won't work very well, so you'll have to use a rational approximation.

If your input numbers are in an arithmetic sequence (some constant step between each consecutive pair), then you can calculate values of the polynomial very quickly by forward differencing.

I'm not sure that this approach is any better than the other one suggested. Depends on what sort of hardware you have, what accuracy you require, etc.

bubba
  • 45,007
  • Can you please share some details on how to use forward differencing to calculate the square roots? – Ecir Hana Jan 26 '13 at 18:13
  • This explains how to it with cubic polynomials: – bubba Jan 27 '13 at 02:48
  • Sorry, I hit the return key by mistake. Here's the link: http://www.drdobbs.com/forward-difference-calculation-of-bezier/184403417 – bubba Jan 27 '13 at 02:49
  • The point is that, once the process gets going, you can calculate the value of a degree $n$ polynomial just by doing $n+1$ additions. If a fairly small value of $n$ gives you the accuracy you need, then this will be very fast. Errors will accumulate at each step, but probably not enough to hurt you. – bubba Jan 27 '13 at 03:04
  • Thanks, but is it possible to calculate "square roots" using FD? – Ecir Hana Jan 27 '13 at 09:41
  • Yes. First construct a polynomial that approximates the square root function over the range that interests you. Once you have a polynomial, you can use FD. – bubba Jan 30 '13 at 07:31
0

As promised, here is another way to compute the square root using rational approximations. The method converges as fast as Newton's method, but is not as well known. I read about this in AMS monthly, but unfortunately I don't remember the issue or authors. If anyone has a pointer to the article that would be great. Let's use the formula $$ \sqrt{n+1}=\frac{n}{\sqrt{n-1}} \sqrt{1-1/n^2} $$ as a starting point. If $\sqrt{n-1}$ is known, we only need to compute $\sqrt{1-1/n^2}$. Let's take $x=-1/n^2$ and let's try to find a fraction $(a+bx)^2/(c+dx)^2$ that is close to $1+x$. It can be checked that $$ \frac{(4+3x)^2}{(4+x)^2}=1+x-\frac{x^3}{16}+\frac{x^4}{32}+\cdots=1+x+\mathcal{O}(x^3) $$ does the trick. Therefore we have $$ \sqrt{1+x}=\frac{4+3x}{4+x}\sqrt{(1+x)\frac{(4+x)^2}{(4+3x)^2}}=\frac{4+3x}{4+x}\sqrt{1+\frac{x^3}{(4+3x)^2}} $$ We can now repeat the whole procedure again by using $\frac{x^3}{(4+3x)^2}$ in place of $x$. Thus we take $x_1=-1/n^2$, $x_2=x_1^3/(4+3x_1)^2$ and we get $$ \sqrt{1-\frac{1}{n^2}}\approx \frac{(4+3x_1)(4+3x_2)}{(4+x_1)(4+3x_2)} $$ and the error is already of order $\mathcal{O}(n^{-18})$, which for $n=10$ should be enough to reconstruct the square-root to double precision. Even the first term $\frac{4+3x_1}{4+x_1}$ is probably pretty good for practical purposes.

Peder
  • 2,252