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.

@MadPhysicist I am sorry, the question I am asking is, is my code wrong or is it actually something with the compiler on Repl.it? — , Sep 05 '18 at 02:23
Are you aware of the fact that no number that has a finite binary representation can be squared to obtain exactly 2? — Mad Physicist, Sep 05 '18 at 02:23
When I run it, it literally just runs without ever producing an output. — , Sep 05 '18 at 02:26
https://stackoverflow.com/questions/3420812/how-do-i-find-if-two-variables-are-approximately-equals , https://stackoverflow.com/questions/14844659/comparing-float-and-double-values-with-delta/ — Mitch Wheat, Sep 05 '18 at 02:34
sqrt is not the best name to use for variable, since it is also a function in math, numpy, and scipy. — ShpielMeister, Sep 05 '18 at 02:36
@ShpielMeister With Python's modules this isn't a big trouble because `sqrt` doesn't interfere with `math.sqrt`. — iBug, Sep 05 '18 at 02:39
yes, unless someone does `from math import *`. in general it's just best to avoid confusion for the reader as well. in your code `i` is the square root, and `sqrt` is the square of putative square root. it is a bit muddled. — ShpielMeister, Sep 05 '18 at 02:44

iBug · Answer 1 · 2018-09-05T02:50:45.690

5

The problem is that using == on floating point numbers is almost always a bad practice and there are v a r i o u s 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)

edited Sep 05 '18 at 02:50

answered Sep 05 '18 at 02:26

iBug

35,554
7
89
134

Might want to mention that precision has to be greater than ~5e-16 if you want the loop to stop. – Mad Physicist Sep 05 '18 at 02:31
I removed the downvote. I generally consider wrong information to warrant a downvote since the purpose is to be helpful to others. I do try to leave meaningful comments though. Thanks for making the fix. – Mad Physicist Sep 05 '18 at 02:36
@MadPhysicist Whoops, I appreciate that though not everyone takes it friendly like me. – iBug Sep 05 '18 at 02:39
I've noticed that. I don't take it friendly when I see mistakes, so I apologize if I was rude. That was not my intent. – Mad Physicist Sep 05 '18 at 03:09
That being said, I've posted an alternative solution. – Mad Physicist Sep 05 '18 at 03:09

Mad Physicist · Answer 2 · 2018-09-05T12:31:52.910

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

I'm not sure your logic is right. You don't necessarily return the integer sqrt even if one exists, and if you set num to 9, you seem to get 1 as the result. — DSM, Sep 05 '18 at 03:33
@DSM. The test after the loop only applies if there was no break. Edited — Mad Physicist, Sep 05 '18 at 03:41
@DSM. I've just removed all the other hard coded references to 2 besides the initial assignment — Mad Physicist, Sep 05 '18 at 03:43
An alternative to checking if you reached machine precision is to iterate 53 times (in the context IEEE double precision). — , Sep 05 '18 at 07:27
@Yves. My method would work with all of the other common float types as well. For example, numpy has float128, float64, float32, etc. Your suggestion, however, would work well if num were an arbitrary-precision type, like a Python int, but for floating point. — Mad Physicist, Sep 05 '18 at 12:19
@MadPhysicist: you needn't tell me. But it is not necessary to recompute the number of bits every time, which is what this test amounts to. (Though the overhead is tiny.) — , Sep 05 '18 at 12:23

pfctgeorge · Answer 3 · 2018-09-05T02:29:52.307

1

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)

edited Sep 05 '18 at 02:29

answered Sep 05 '18 at 02:24

pfctgeorge

698
3
9

ShpielMeister · Answer 4 · 2018-09-05T03:58:34.840

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)

score 0 · Answer 5 · 2018-09-05T07:23:09.933

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.

Binary search for square root 2

5 Answers5