Square root with babylonian or Herons method?

Question

Im implementing in java the Babylonian/Heron`s method to get the square root of a number, based on Wikipedia info

Currently i have:

public static void main(String[] args) {
    System.out.println(Sqrt.sqrt1(8));
    //System.out.println(Sqrt.sqrt2(8));   //Infinite loop
    System.out.println(Sqrt.sqrt3(8));
    System.out.println(Sqrt.sqrt4(8));
}

static float sqrt1(float x) {
    float b = 0, h = x;

    while (b != h) {
        b = (h + b) / 2;
        h = x / b;
    }
    return b;
}

static double sqrt2(double x) {
    double b = x, h = 0;

    while (b != h) {
        b = (h + b) / 2;
        h = x / b;
    }
    return b;
}

static double sqrt3(double x) {
    double b = x, h = 0;

    while (Math.abs(b - h) > 0.000000000001) {
        b = (h + b) / 2;
        h = x / b;
    }
    return b;
}

static double sqrt4(double x) {
    double r = x, t = 0;

    while (t != r) {
        t = r;
        r = (x / r + r) / 2;
    }
    return r;
}

The output would be:

2.828427

2.82842712474619

2.82842712474619

The sqrt2 method will loop forever, but this only happens with doubles because the sqrt1 method works fine with floats. I dont know why is this. So the sqrt3 method looks like the way to go if i want to work with doubles.

Im a bit confused about which methods im implementing, ¿Is the Babylonian method the same as the Herons method?.

From what i understand, the Babylonian method is based upon the fact that a side of a square is the square root of area of the square (x). So you can start with a rectangle with b.h dimensions, get the average of the two sides (b=b+h/2) and then consider this result as the side of a smaller rectangle and of course get the other side (h=x/b). The rectangle will start aproximating into the desired square. This is what i have done in sqrt1, sqrt2 and sqrt3 methods:

while (b != h) {
    b = (h + b) / 2;
    h = x / b;
}

In the other side, the Wikipedia link says the Babylonian/Herons method is the same and describe it as:

"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"

You can see this implementation in the sqrt4 method:

    while (t != r) {  
        t = r;  
        r = (x/r + r) / 2;  
    }  

From what i can see, this two methods are not the same, but similar. Please correct me if im wrong.

Something interesting is that i cant do: while (b != h) { when b and h are doubles like showed in the sqrt2 method because it will loop forever. Instead i used while (Math.abs(b - h) > 0.000000000001) {

However, i can do: while (t != r) { with b and h doubles as showed in sqrt4.

I would appreciate somebody explain this behaivour.

-----------EDIT:------------------------------------------------------------

As suggested, both algorithms are the same but implemented diferent. However the following suggested code loops for ever just like sqrt2 with x=8:

static double sqrt5(double x) {
    double b = x;

    while (b != x/b) {
        b = (x / b + b) / 2;
    }
    return b;
}

So.. what is the difference in sqrt4 with the other implementations? Why it doesnt loop forever like the others?

Answer 1

Update: sqrt4 has better chances of eventually exiting the loop than sqrt5 because its stopping condition compares the approximation from one iteration with one from the iteration before it.

The computation tends to reduce the error, so that eventually the value computed for b is so close to the exact square root that x/b differs from b only in the last bit. At this point the value computed for (x / b + b) / 2 using the finite precision available will equal b, and the iteration stops.

For example, if we are calculating sqrt(2) and have reached the approximation b = 1.414213562373095, we have:

>>> b
1.414213562373095
>>> 2/b                     # Close but not quite equal to b, 
1.4142135623730951          # iteration in sqrt5 continues
>>> (2/b + b)/2        
1.414213562373095           # Exactly equal to b, sqrt4 stops

As you can see, once b has reached 1.414213562373095 its value will not be changed by the calculation in the loop. Because b and 2/b still differ in the last digit sqrt5 will never exit.

---

The Babylonian and Heron's methods are the same algorithm, and it coincides with the Newton-Rhapson method for solving x²=a. The difference between the various implementations you have is the stopping condition. For example:

while (b != h) {
    b = (h + b) / 2;
    h = x / b;
}

is the same as:

while (b != x/b) {
    b = (x / b + b) / 2;
}

Of course, b != x/b is not a very good stopping condition because b may never become the mathematically exact square root of x: if b is not the exact square root, x/b does not equal b. For example in Python:

>>> sqrt(2)
1.4142135623730951
>>> 2/sqrt(2)
1.4142135623730950

You should almost always use a bound on the relative difference as the stopping condition, for example

eps = 1e-6;
while (abs(b-h)/b > eps) { ...

Answer 2

Because of rounding errors, you should never rely on equality of doubles or floats. You need another method for determining when your approximation will be good enough.

A good start is to print out the approximation in every cycle so you see how it goes. It may be that the difference between two consecutive approximations gets smaller and smaller and then, suddenly, becomes chaotic. Then you know that you did go too far, and go back to the last approcximation where the difference to the previous was smaller than the previous difference.

Be sure to test cases that would yield a value that can never be represented in a double, like sqrt(2). Be sure to test cases that should yield a number that is representable. Finally, make sure your function will yield an integral value when the square root of a square is taken.

For the differences with float and double: as floats have less precision, you reach a point where the difference is too small and due to rounding errors results in 0. But this is mere luck, and could be different with different inputs. However, the less bits you have in the mantissa, the more likely it is.

Answer 3

As others have written, the Babylonian method and Heron's method are the same thing. The Babylonians invented it, but Heron was the first to write it down, several hundred years later, so he gets the credit. Sometimes life isn't fair.

I also agree with the comments about precision and rounding errors. But I have a different solution to propose. Normalize the number x for which the square root is being sought onto the range 1 <= x < 4 by continually multiplying or dividing by 4 until it is in range. Then "unroll the loop" that performs the calculations; since x is in a known range, the number of loops can be calculated in advance, without the need to make a calculation to determine when to terminate the loop. Finally, multiply or divide by 2 as many times as you divided or multiplied by 4 initially. Here's the code:

function sqrt(n)
    if n < 1
        return sqrt(n * 4) / 2
    if 4 <= n
        return sqrt(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'll leave it to you to translate the pseudocode above to Java. Five unrollings of the loop are sufficient for doubles. Since each time through the loop doubles the number of digits of accuracy, four unrollings of the loop are sufficient for floats.

I discuss Heron's method at my blog.