Mathematics behind Babylonian Square Root method [

2020-04-20 09:50发布

问题:


Want to improve this question? Update the question so it's on-topic for Stack Overflow.

Closed 6 years ago.

I read the method to calculate the square root of any number and the algorithm is as follows:

double findSquareRoot(int n) {
    double x = n;
    double y = 1;
    double e = 0.00001;
    while(x-y >= e) {
        x = (x+y)/2;
        y = n/x;
    }
    return x;
}

My question regarding this method are

  1. How it calculates the square root? I didn't understand the mathematics behind this. How x=(x+y)/2 and y=n/x converges to square root of n. Explain this mathematics.

  2. What is the complexity of this algorithm?

回答1:

It is easy to see if you do some runs and print the successive values of x and y. For example for 100:

50.5 1.9801980198019802
26.24009900990099 3.8109612300726345
15.025530119986813 6.655339226067038
10.840434673026925 9.224722348894286
10.032578510960604 9.96752728032478
10.000052895642693 9.999947104637101
10.000000000139897 9.999999999860103

See, the trick is that if x is not the square root of n, then it is above or below the real root, and n/x is always on the other side. So if you calculate the midpoint of x and n/x it will be somewhat nearer to the real root.

And about the complexity, it is actually unbounded, because the real root will never reached. That's why you have the e parameter.



回答2:

This is a typical application of Newton's method for calculating the square root of n. You're calculating the limit of the sequence:

x_0 = n
x_{i+1} = (x_i + n / x_i) / 2

Your variable x is the current term x_i and your variable y is n / x_i.

To understand why you have to calculate this limit, you need to think of the function:

f(x) = x^2 - n

You want to find the root of this function. Its derivative is

f'(x) = 2 * x

and Newton's method gives you the formula:

x_{i+1} = x_i - f(x_i) / f'(x_1) = ... = (x_i + n / x_i) / 2

For completeness, I'm copying here the rationale from @rodrigo's answer, combined with my comment to it. This is helpful if you want to forget about Newton's method and try to understand this algorithm alone.

The trick is that if x is not the square root of n, then it is an approximation which lies either above or below the real root, and y = n/x is always on the other side. So if you calculate the midpoint of (x+y)/2, it will be nearer to the real root than the worst of these two approximations (x or y). When x and y are close enough, you're done.

This will also help you find the complexity of the algorithm. Say that d is the distance of the worst of the two approximations to the real root r. Then the distance between the midpoint (x+y)/2 and r is at most d/2 (it will help you if you draw a line to visualize this). This means that, with each iteration, the distance is halved. Therefore, the worst-case complexity is logarithmic w.r.t. to the distance of the initial approximation and the precision that is sought. For the given program, it is

log(|n-sqrt(n)|/epsilon)


回答3:

I think all information can be found in wikipedia.

The basic idea is that if x is an overestimate to the square root of a non-negative real number S then S/x, will be an underestimate and so the average of these two numbers may reasonably be expected to provide a better approximation.

With each iteration this algorithm doubles correct digits in answer, so complexity is linear to desired accuracy's logarithm.

Why does it work? As stated here, if you will do infinite iterations you'll get some value, let's name it L. L has to satisfy equasion L = (L + N/L)/2 (as in algorithm), so L = sqrt(N). If you're worried about convergence, you may calculate squared relative errors for each iteration (Ek is error, Ak is computed value):

Ek = (Ak/sqrt(N) - 1)²
if:
Ak = (Ak-1 + N/Ak-1)/2 and Ak = sqrt(N)(sqrt(Ek) + 1)
you may derive recurrence relation for Ek:
Ek = Ek-1²/[4(sqrt(Ek-1) + 1)²]
and limit of it is 0, so limit of A1,A2... sequence is sqrt(N).



回答4:

The mathematical explanation is that, over a small range, the arithmetic mean is a reasonable approximation to the geometric mean, which is used to calculate the square root. As the iterations get closer to the true square root, the difference between the arithmetic mean and the geometric mean vanishes, and the approximation gets very close. Here is my favorite version of Heron's algorithm, which first normalizes the input n over the range 1 ≤ n < 4, then unrolls the loop for a fixed number of iterations that is guaranteed to converge.

def root(n):
    if n < 1: return root(n*4) / 2
    if 4 <= n: return root(n/4) * 2
    x = (n+1) / 2
    x = (x + n/x) / 2
    x = (x + n/x) / 2
    x = (x + n/x) / 2
    x = (x + n/x) / 2
    x = (x + n/x) / 2
    return x

I discuss several programs to calculate the square root at my blog.