sqrt function weird behaviour in Python

Question

In Python, I wrote a custom code sqrt(x, delta) to calculate the square root of a given number with a delta-close approximation. It uses a while loop and a binary-search-like algorithm.
The code:

from __future__ import division

def sqrt(x, delta):
    start = 0
    end = x
    while (end-start) > delta:
        middle = (start + end) / 2
        middle_2 = middle * middle
        if middle_2 < x:
            start = middle
            print "too low"
        elif middle_2 > x:
            end = middle
            print "too high"
        else:
            return middle
    result = (start + end) / 2
    return result

It basicly works and it is quite fast, but there are cases when it gets into an infinite while-loop.
Examples:

sqrt(1e27, 1/1024) => works fine (returns 'too low's and 'too high's, then returns correct result)
sqrt(1e28, 1/1024) => works fine
sqrt(1e29, 1/1024) => never-ending loop, it keeps printing 'too low' forever
sqrt(1e30, 1/1024) => works fine
sqrt(1e31, 1/1024) => 'too low' forever
sqrt(1e32, 1/1024) => works fine
sqrt(1e33, 1/1024) => works fine (also surprising after the problem with 1e29 and 1e31)
sqrt(1e34, 1/1024) => works fine
sqrt(1e35, 1/1024) => 'too low' forever
sqrt(1e36, 1/1024) => works fine
sqrt(1e37, 1/1024) => 'too high' forever (too high this time!)
sqrt(1e38, 1/1024) => works fine
sqrt(1e39, 1/1024) => works fine (surprising again..)
... 1e40-1e45 ... they all work fine
sqrt(1e46, 1/1024) => 'too low' forever (surprisingly it occurs now with 1e'even number')
...
sqrt(1e200, 1/1024) => works fine
sqrt(1e201, 1/1024) => works fine
...
sqrt(1e299, 1/1024) => 'too low' forever
sqrt(1e300, 1/1024) => 'too high' forever
...
sqrt(1e304, 1/1024) => 'too high' forever
sqrt(1e305, 1/1024) => works fine
... 305-308 ... they allwork fine
sqrt(1e309, 1/1024) => inf (reached some 'infinite' limit?)

I first thought it was with numbers above a limit, like 1e20.. But then it worked with them as well. Also, I had this thought that is was about 1e'odd' or 1e'even' numbers, but as we can see in the examples, it wasn't that. I also tried it using different delta's instead 1/1024, but they showed similar behavior.

I would appreciate any explanation that tells what is going on behind the scenes that causes this behavior.

Answer 1

float can only represent a finite set of numbers. Your code ends up in a situation where start and end are two consecutive such numbers. As a result, (start + end) / 2 has to be either rounded down to start or rounded up to end.

If it gets rounded down, middle_2 < x. Now, if end - start > delta, you've got a "too low" infinite loop.

If it gets rounded up, and if end - start > delta, you've got a "too high" infinite loop.

You should probably redefine delta as a relative error rather than an absolute one.