2
$\begingroup$

I implemented Newton's method to find a square root. I then realized that the Babylonian method (which is just another way to look at Newton's method) looks a lot like a binary search. I then implemented a binary search method to compute the SQRT of a number.

I can't find online any examples of binary search as a valid algorithm for Square root (the wikipedia page doesn't mention anything, and I haven't found an example here on math.stackexchange...)

I thought the binary search should be just as fast, but at lower values of N (i.e. 2) it is 10x slower... For bigger values of N (i.e. N=234567919191) it is 3x slower. For huge values of N (i.e. N=10^154) it is 1.16x slower. This makes me wonder if there is a N big enough where Binary Search is equal to Newton/Babylonian? (could it ever be better?)

I don't understand why though. If I use the "guess" or "seed" to be 1, we have the following:

  • Binary Search: mean(low,high), with starting low=1, high=n, then adjust our low/high boundaries to [1,mid] or [mid,N] as needed, and iterate over until the relative (or absolute) error of our mid is within our parameters.
  • Babylonian Search: mean(guess, N/guess) and iterate over. But if I guess 1, then we have: mean(1, N) which is the same exact starting point as the binary search! But then the second call would be: mean((N+1)/2, ((N+1)/2)/N).

If we look at N=100, the first few iterations of each method:

Binary Search:

low high mean (of high-low)
0 100 50.0
0 50.0 25.0
0 25.0 12.5
0 12.5 6.25
6.25 12.5 3.125
9.375 12.5 1.5625
9.375 10.9375 0.78125
9.375 10.15625 0.390625

Babylonian Method:

low high mean
1 100.0 50.5
1.9801980198019802 50.5 26.24009900990099
3.8109612300726345 26.24009900990099 15.025530119986813
6.655339226067038 15.025530119986813 10.840434673026925
9.224722348894286 10.840434673026925 10.032578510960604
9.96752728032478 10.032578510960604 10.000052895642693
9.999947104637101 10.000052895642693 10.000000000139897
9.999999999860103 10.000000000139897 10.0

So it looks like the low/high boundaries of the Babylonian Search converge faster than just a pure binary search. I just don't know how to show this mathematically... Maybe someone can help? This is for my own curiosity, nothing more.

I guess I could adjust the binary search so the boundary that is not moving (low or high) gets adjusted as well, so it converges faster? But then isn't it just the Babylonian method again? :D

I have the code on github: https://github.com/CobraCoral/recursive_square/blob/master/compute_sqrt.py I also pasted it below for your convenience, thanks!


Newton's method:

def sqrt(N, X=None, iterations=0):
    """
    Newton's method
    """
    if N < 0:
        result = sqrt(-N, X, iterations)
        return (complex(result[0], 1), result[1])
    if N == 0: return (0, iterations)
    if N == 1: return (1, iterations)
    if not X: X = N/2 # seeding our first X ... it can be any number [0,N)
    while abs(1 - (X**2 / N)) > sqrt.error:  # relative error
        X = X - ((X**2 - N) / (2 * X))
        iterations += 1
    return (X, iterations)
sqrt.error = 0.00000000000001

Babylonian SQRT method:

def babylonian_sqrt(N, X=None, iterations=0):
    """
    Babylonian method:
        Xi: Some guess between [0, N)
        N: Number we want to find the square root of
        Formula: (Xi + (N / Xi)) / 2  -> This is also the arithmetic mean of the guess, [Xi, N/Xi]
    """
    if N < 0:
        result = babylonian_sqrt(-N, iterations)
        return (complex(result[0], 1), result[1])
    if N == 0: return (0, iterations)
    if N == 1: return (1, iterations)
    if not X: X = N/2 # seeding our first X ... it can be any number [0,N)
    while abs(1 - (X**2 / N)) > babylonian_sqrt.error:  # relative error
        X = (X + (N / X)) / 2
        iterations += 1
    return (X, iterations)
babylonian_sqrt.error = 0.00000000000001

Binary Search method:

def sqrt_binary_search(N, low=0, high=None, iterations=0):
    """
    Binary search method:
    """
    if N < 0:
        result = sqrt_binary_search(-N, low, high, iterations)
        return (complex(result[0], 1), result[1])
    if N == 0: return (0, iterations)
    if N == 1: return (1, iterations)
    high = N if not high else high
    mid = 0
    while abs(1 - ((low + mid)**2/N)) > sqrt_binary_search.error:  # relative error
        mid = (high - low) / 2
        newX_square = (low + mid) ** 2
        if newX_square > N:
            high = high - mid
        else:
            low = low + mid
        iterations += 1
    return (low + mid, iterations)
sqrt_binary_search.error = 0.00000000000001
$\endgroup$

1 Answer 1

1
$\begingroup$

In binary search, the length of the interval is multiplied by $1/2$ on each iteration. But in Newton's method (when you're close to the true root, and the derivative of the function at the root is nonzero) the error $\epsilon_n$ at iteration $n$ satisfies a bound of the form $\epsilon_{n+1} \le c \epsilon_n^2$. When $\epsilon_n$ is small this converges much faster than binary search.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .