Binary search for square root 2

Question

I have an issue where my binary search algorithm to find the square root of 2 appears to be in an infinite loop and runs forever:

num = 2
low = 1
high = num
i = 0
while((i**2) != 2): #was while(low<high): but wasnt working with it either
 i = low + (high - low) / 2;

 sqrt = i * i

 if (sqrt == num):
     print(i)
 elif(sqrt < num):
     low = i
 else:
     high = i
print(sqrt)    

testRoot = 2 ** (.5)
print(testRoot)

I am not sure if there is an issue with my while loop or not. I assumed it would be a pretty straight forward binary search algorithm with a slight modification to accommodate the square root aspect.

When running my code, I can't seem to get it to produce an output whatsoever. I am unsure if there is a real issue with the code or the compiler as I thought I followed the algorithm pretty closely to what I have in the past.

Answer 1

As mentioned in my original comment and all the answers, the square root of 2 is irrational. The square root of every integer that's not a perfect square is irrational for that matter, so 2 isn't special in that regard. What matters is that x**2 == 2 will never be true for any x of finite precision (since finite precision is another way of saying a number is rational).

The other answers suggest searching until you reach some fixed, pre-determined accuracy. That works well, especially if you know the binary order of magnitude of the answer ahead of time, since then you can set the accuracy of the result to be in the last digit.

I'd like to suggest a more natural approach. You can check if your center value exactly equals one of the bounds. That would mean that half the difference between the bounds represents less than one digit of precision in your current guess. Your phrasing of the center is already correct: i = low + (high - low) / 2 can be compared to low and high using ==, while i = (low + high) / 2 may not. This is because the precision of high - low is greater than or equal to the precision of either bound, while low + high is likely to lose some digits.

So here is what I would recommend:

num = 2
low = 1
high = num
guess = low + (high - low) / 2
count = 0
while guess != low and guess != high:
    sqr = guess * guess

    if sqr == num:
        break
    elif(sqr < num):
        low = guess
    else:
        high = guess

    guess = low + (high - low) / 2
    count += 1
else:
    if abs(low * low - num) < abs(high * high - num):
        guess = low
    else:
        guess = high

print(count, ':', sqr)
print(num ** (.5), '~=', guess)

I've added count for verification. The result is obtained in 52 iterations, accurate within 1 digit of precision:

52 : 2.0000000000000004
1.4142135623730951 ~= 1.4142135623730951 

The final check against the bounds (the else clause of the while) ensures that you get the closest one to the desired result, regardless of which one you hit first.

The convergence is sane: a 64-bit floating point number in IEEE-754 format has 53 bits in the mantissa, so it makes sense that you would have to halve your search space exactly that many times to get your result (the first time is outside the loop).

Here's the snippet I used for testing: https://ideone.com/FU9r82

Answer 2

As said in several other posts, the comparison i * i == 2 cannot work.

A simple solution to design the stopping criterion is to tell how many bits of accuracy you want. Indeed, in a dichotomic search the number of exact bits increases by one on very iteration.

Thus for full double-precision accuracy, iterate 53 times (more is useless). Also note that testing for equality inside the loop is counter-productive.

num= 2
low= 1
hig= 2
for i in range(53):
    mid= 0.5 * (low + hig)
    if mid * mid < num:
        low= mid
    else:
        hig= mid

print(mid)

53 iterations are appropriate here because the initial estimates are accurate for the first bit (1≤√2<2). For less accurate initial estimates, add a few iterations.

Answer 3

It's often helpful to do some exploration.

$ python
Python 3.6.6 
>>> import math

>>> import numpy
>>> import scipy
>>> import numpy
>>> math.sqrt(2) ** 2
 2.0000000000000004
>>> numpy.sqrt(2) ** 2
 2.0000000000000004
>>> scipy.sqrt(2) ** 2
 2.0000000000000004
>>> (2.0**(0.5))**2
 2.0000000000000004
>>> x =  math.sqrt(2) ** 2
>>> math.sqrt(x)
 1.4142135623730951
>>> math.sqrt(2)
 1.4142135623730951
>>> x*x
 4.000000000000002
>>> x**2
 4.000000000000002
>>> 1.414213562373095**2
 1.9999999999999996
>>> 1.41421356237309505**2
 2.0000000000000004
>>> 1.41421356237309505
 1.4142135623730951
>>> 1.41421356237309504
 1.4142135623730951
>>> 1.41421356237309504**2
 2.0000000000000004
>>> 1.41421356237309503**2
 1.9999999999999996
>>> 1.41421356237309503 * 1.41421356237309504
 2.0
>>> 1.41421356237309504 - 1.41421356237309503
 2.220446049250313e-16

A little bit of elbow grease can be educational. (and confusing!)

let's take a look at the errors

>>> s =set()
>>> exact = 0
>>> exact_over = 0
>>> exact_under = 0

>>>for i in range(100):
...     differ = (i - (i**0.5)**2)
...     s |= {differ}
...     exact += 1 if differ == 0 else 0
...     exact_over  += 1 if differ > 0  else 0
...     exact_under += 1 if differ < 0  else 0

>>> sorted(list(s)) 
[-1.4210854715202004e-14,
 -7.105427357601002e-15,
 -3.552713678800501e-15,
 -1.7763568394002505e-15,
 -8.881784197001252e-16,
 -4.440892098500626e-16,
 0.0,
 4.440892098500626e-16,
 8.881784197001252e-16,
 1.7763568394002505e-15,
 3.552713678800501e-15,
 7.105427357601002e-15,
 1.4210854715202004e-14]

>>> exact_under, exact, exact_over
(26, 49, 25)

Answer 4

The problem is that using == on floating point numbers is almost always a bad practice and there are various questions about it. You should replace the comparison with abs(a - b) < precision. Also read the comments under this question, very useful.

My revised code looks like this (Python 3), replace 1e-6 with a smaller number if you want more precision. But note that it's not sensible to be "too precise", and a precision of 1.0e-15 or bigger is recommended if you want the loop to stop, because floating point number itself has precision limits.

num = 2
low = 1
high = num
i = 0
while abs((i**2) - num) > 1e-6:  # Note
    i = low + (high - low) / 2
    sqrt = i * i

    if abs(sqrt - num) < 1e-6:  # Note
        print(i)
    elif(sqrt < num):
        low = i
    else:
        high = i
print(sqrt)

testRoot = 2 ** (.5)
print(testRoot)

Answer 5

Equal between two floats is a very rigid condition because square root of 2 has a infinite number of digit after decimal point. Try this while condition:

while (abs((i ** 2) - 2) > 1e-8)