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

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.)

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).