Inverse function for monotonically increasing function, OverflowError for log10()

Question

For an assignment, we were asked to create a function which returns an inverse function. The basic problem was to create a square root function from a square function. I came up with a solution using binary search and another solution using Newton's method. My solution seems to work fine for cube-root and square-root but not for log10. Here are my solutions:

#Binary Search
def inverse1(f, delta=1e-8):
    """Given a function y = f(x) that is a monotonically increasing function on
    non-negative numbers, return the function x = f_1(y) that is an approximate
    inverse, picking the closest value to the inverse, within delta."""
    def f_1(y):
        low, high = 0, float(y)
        last, mid = 0, high/2
        while abs(mid-last) > delta:
            if f(mid) < y:
                low = mid
            else:
                high = mid
            last, mid = mid, (low + high)/2
        return mid
    return f_1

#Newton's Method
def inverse(f, delta=1e-5):
    """Given a function y = f(x) that is a monotonically increasing function on
    non-negative numbers, return the function x = f_1(y) that is an approximate
    inverse, picking the closest value to the inverse, within delta."""
    def derivative(func): return lambda y: (func(y+delta) - func(y)) / delta
    def root(y): return lambda x: f(x) - y
    def newton(y, iters=15):
        guess = float(y)/2
        rootfunc = root(y)
        derifunc = derivative(rootfunc)
        for _ in range(iters):
            guess = guess - (rootfunc(guess)/derifunc(guess))
        return guess
    return newton

Regardless which method is used, when I get to the input n = 10000 for log10() in the professor's test function, I get this error: (EXCEPTION: when my newton's method function is used, log10() is way off, whereas this binary-search method is relatively accurate until the input threshold is reached, either way, both solutions throw this error when n = 10000)

   2: sqrt =     1.4142136 (    1.4142136 actual); 0.0000 diff; ok
   2: log =     0.3010300 (    0.3010300 actual); 0.0000 diff; ok
   2: cbrt =     1.2599211 (    1.2599210 actual); 0.0000 diff; ok
   4: sqrt =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
   4: log =     0.6020600 (    0.6020600 actual); 0.0000 diff; ok
   4: cbrt =     1.5874011 (    1.5874011 actual); 0.0000 diff; ok
   6: sqrt =     2.4494897 (    2.4494897 actual); 0.0000 diff; ok
   6: log =     0.7781513 (    0.7781513 actual); 0.0000 diff; ok
   6: cbrt =     1.8171206 (    1.8171206 actual); 0.0000 diff; ok
   8: sqrt =     2.8284271 (    2.8284271 actual); 0.0000 diff; ok
   8: log =     0.9030900 (    0.9030900 actual); 0.0000 diff; ok
   8: cbrt =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
  10: sqrt =     3.1622777 (    3.1622777 actual); 0.0000 diff; ok
  10: log =     1.0000000 (    1.0000000 actual); 0.0000 diff; ok
  10: cbrt =     2.1544347 (    2.1544347 actual); 0.0000 diff; ok
  99: sqrt =     9.9498744 (    9.9498744 actual); 0.0000 diff; ok
  99: log =     1.9956352 (    1.9956352 actual); 0.0000 diff; ok
  99: cbrt =     4.6260650 (    4.6260650 actual); 0.0000 diff; ok
 100: sqrt =    10.0000000 (   10.0000000 actual); 0.0000 diff; ok
 100: log =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
 100: cbrt =     4.6415888 (    4.6415888 actual); 0.0000 diff; ok
 101: sqrt =    10.0498756 (   10.0498756 actual); 0.0000 diff; ok
 101: log =     2.0043214 (    2.0043214 actual); 0.0000 diff; ok
 101: cbrt =     4.6570095 (    4.6570095 actual); 0.0000 diff; ok
1000: sqrt =    31.6227766 (   31.6227766 actual); 0.0000 diff; ok
Traceback (most recent call last):
  File "/CS212/Unit3HW.py", line 296, in <module>
    print test()
  File "/CS212/Unit3HW.py", line 286, in test
    test1(n, 'log', log10(n), math.log10(n))
  File "/CS212/Unit3HW.py", line 237, in f_1
    if f(mid) < y:
  File "/CS212/Unit3HW.py", line 270, in power10
    def power10(x): return 10**x
OverflowError: (34, 'Result too large')

Here is the test function:

def test():
    import math
    nums = [2,4,6,8,10,99,100,101,1000,10000, 20000, 40000, 100000000]
    for n in nums:
        test1(n, 'sqrt', sqrt(n), math.sqrt(n))
        test1(n, 'log', log10(n), math.log10(n))
        test1(n, 'cbrt', cbrt(n), n**(1./3.))


def test1(n, name, value, expected):
    diff = abs(value - expected)
    print '%6g: %s = %13.7f (%13.7f actual); %.4f diff; %s' %(
        n, name, value, expected, diff,
        ('ok' if diff < .002 else '**** BAD ****'))

These is how the test is setup:

#Using inverse() or inverse1() depending on desired method
def power10(x): return 10**x
def square(x): return x*x
log10 = inverse(power10)
def cube(x): return x*x*x
sqrt = inverse(square)
cbrt = inverse(cube)
print test()

The other solutions posted seem to have no problems running the full set of test inputs (I have tried not to look at the posted solutions). Any insight into this error?

---

It seems as though the consensus is the size of the number, however, my professor's code seems to work just fine for all cases:

#Prof's code:
def inverse2(f, delta=1/1024.):
    def f_1(y):
        lo, hi = find_bounds(f, y)
        return binary_search(f, y, lo, hi, delta)
    return f_1

def find_bounds(f, y):
    x = 1
    while f(x) < y:
        x = x * 2
    lo = 0 if (x ==1) else x/2
    return lo, x

def binary_search(f, y, lo, hi, delta):
    while lo <= hi:
        x = (lo + hi) / 2
        if f(x) < y:
            lo = x + delta
        elif f(x) > y:
            hi = x - delta
        else:
            return x;
    return hi if (f(hi) - y < y - f(lo)) else lo

log10 = inverse2(power10)
sqrt = inverse2(square)
cbrt = inverse2(cube)

print test() 

RESULTS:

     2: sqrt =     1.4134903 (    1.4142136 actual); 0.0007 diff; ok
     2: log =     0.3000984 (    0.3010300 actual); 0.0009 diff; ok
     2: cbrt =     1.2590427 (    1.2599210 actual); 0.0009 diff; ok
     4: sqrt =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
     4: log =     0.6011734 (    0.6020600 actual); 0.0009 diff; ok
     4: cbrt =     1.5865107 (    1.5874011 actual); 0.0009 diff; ok
     6: sqrt =     2.4486818 (    2.4494897 actual); 0.0008 diff; ok
     6: log =     0.7790794 (    0.7781513 actual); 0.0009 diff; ok
     6: cbrt =     1.8162270 (    1.8171206 actual); 0.0009 diff; ok
     8: sqrt =     2.8289337 (    2.8284271 actual); 0.0005 diff; ok
     8: log =     0.9022484 (    0.9030900 actual); 0.0008 diff; ok
     8: cbrt =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
    10: sqrt =     3.1632442 (    3.1622777 actual); 0.0010 diff; ok
    10: log =     1.0009756 (    1.0000000 actual); 0.0010 diff; ok
    10: cbrt =     2.1534719 (    2.1544347 actual); 0.0010 diff; ok
    99: sqrt =     9.9506714 (    9.9498744 actual); 0.0008 diff; ok
    99: log =     1.9951124 (    1.9956352 actual); 0.0005 diff; ok
    99: cbrt =     4.6253061 (    4.6260650 actual); 0.0008 diff; ok
   100: sqrt =    10.0004883 (   10.0000000 actual); 0.0005 diff; ok
   100: log =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
   100: cbrt =     4.6409388 (    4.6415888 actual); 0.0007 diff; ok
   101: sqrt =    10.0493288 (   10.0498756 actual); 0.0005 diff; ok
   101: log =     2.0048876 (    2.0043214 actual); 0.0006 diff; ok
   101: cbrt =     4.6575475 (    4.6570095 actual); 0.0005 diff; ok
  1000: sqrt =    31.6220242 (   31.6227766 actual); 0.0008 diff; ok
  1000: log =     3.0000000 (    3.0000000 actual); 0.0000 diff; ok
  1000: cbrt =    10.0004883 (   10.0000000 actual); 0.0005 diff; ok
 10000: sqrt =    99.9991455 (  100.0000000 actual); 0.0009 diff; ok
 10000: log =     4.0009756 (    4.0000000 actual); 0.0010 diff; ok
 10000: cbrt =    21.5436456 (   21.5443469 actual); 0.0007 diff; ok
 20000: sqrt =   141.4220798 (  141.4213562 actual); 0.0007 diff; ok
 20000: log =     4.3019052 (    4.3010300 actual); 0.0009 diff; ok
 20000: cbrt =    27.1449150 (   27.1441762 actual); 0.0007 diff; ok
 40000: sqrt =   199.9991455 (  200.0000000 actual); 0.0009 diff; ok
 40000: log =     4.6028333 (    4.6020600 actual); 0.0008 diff; ok
 40000: cbrt =    34.2003296 (   34.1995189 actual); 0.0008 diff; ok
 1e+08: sqrt =  9999.9994545 (10000.0000000 actual); 0.0005 diff; ok
 1e+08: log =     8.0009761 (    8.0000000 actual); 0.0010 diff; ok
 1e+08: cbrt =   464.1597912 (  464.1588834 actual); 0.0009 diff; ok
None

Crosier

Yes, 10**1000.0 is too large. – Shantell

@KennyTM the other students' solutions seem to pass the tests. I feel like my solution is missing something. – Crosier

Perhaps you shouldn't start from y=1000. Try y=1 for all cases. – Shantell

@KennyTM the test cases all start from 2 and increment upwards, and I have tested them for y=1 and they seem to work well for those basic cases. – Crosier

Here is a hint why your professors code is running better. Take a look at the purpose of his "find bounds" method and see why you are getting an overflow – Visually

For successive approximation algorithms like these, a good-enough initial guess is often half the battle. Binary search is like a Aesop's tortoise: slow and dependable if lo and hi bracket the solution. Newton's method is more like a hare. It may quickly jump to the solution (good initial guess), but it may jump off toward infinity (bad initial guess). – Yoga

Answer 1

This is actually a problem in your understanding of the math instead of the program. The algorithm is fine, but the supplied initial condition is not.

You define inverse(f, delta) like this:

def inverse(f, delta=1e-5):
    ...
    def newton(y, iters=15):
        guess = float(y)/2
        ...
    return newton

so you guess the result of 1000 = 10^x is 500.0, but surely 10⁵⁰⁰ is too large. The initial guess should chosen to be in a valid for f, not chosen for the inverse of f.

I suggested that you initialize with a guess of 1, i.e. replace that line with

guess = 1

and it should work fine.

---

BTW, your binary search's initial condition is also wrong, because you assume the solution is between 0 and y:

low, high = 0, float(y)

this is true for your test cases, but it's easy to construct counter examples e.g. log₁₀ 0.1 (= -1), √0.36 (= 0.6), etc. (Your professor's find_bounds method does solve the √0.36 problem, but still won't handle the log₁₀ 0.1 problem.)

Shantell

his prof's find bounds method attempts to solve your second point. – Visually

I guess I automatically assumed that, this being a general inverse function for any mon increasing func, splitting the range down the middle would be a good starting point... should've given it more thought. So staring at 1 would be a better starting point in general or just for log10? – Crosier

Also, any insight as to why my newton's method is way off for all log10 tests (example above was with my bin-search method)? Or does that have to do with the starting guess issue with newton? – Crosier

@Crosier I guess it's better ask on math.stackexchange.com for detail about these numerical root-finding methods. – Shantell

Answer 2

I traced your error but it basically comes down to the fact that 10**10000000 causes an overflow in python. a quick check using the math library

math.pow(10,10000000)

Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    math.pow(10,10000000)
OverflowError: math range error

I did a little research for you and found this

Handling big numbers in code

You need to either re-evaluate why you need to calculate such a large number (and change your code accordingly :: suggested) or start looking for some even larger number handling solutions.

You could edit your inverse function to check if certain inputs will cause it to fail (try statement) which could also solve some issues with zero division if the functions arent monotonically increasing, and avoid those regions OR

you could mirror a number of points in the "interesting" area about y=x and use an interpolation scheme through those points to create the "inverse" function (hermite's, taylor series, etc).

Visually

I wasn't aware of the size restrictions, but it doesn't seem to effect the other solutions posted in the class forum, those pass all tests. Also, any insight into why my Newton's method can't handle log10 at all (answers are wayyyy off)? Thanks – Crosier

The math module essentially provides access to C's fast <math.h> and hence will raise OverflowError exceptions for values where the built-in pow() or the ** operator wouldn't. That said, your larger point is spot-on, as calculating such huge numbers would be prohibitively slow even if it didn't raise an exception. – Croft

Here is a hint why your professors code is running better. Take a look at the purpose of his "find bounds" method and see why you are getting an overflow. You are probably not even getting through one calculation – Visually

Answer 3

If you're using Python 2.x, int and long are distinct types, and the OverflowError can only result for ints (q.v., Built-in Exceptions). Try explicitly using longs instead (either by using the long() built-in function on your integral values, or appending L to the numeric literals).

Edit: Obviously, as Paul Seeb and KennyTM have pointed out in their superior answers, this is no remedy for an algorithmic flaw.

Croft