python - 单调递增函数的反函数，log10() 的 OverflowError

Question

对于一项作业，我们被要求创建一个返回反函数的函数。基本问题是从平方函数创建平方根函数。我想出了一个使用二进制搜索的解决方案和另一个使用牛顿法的解决方案。我的解决方案似乎适用于立方根和平方根，但不适用于 log10。这是我的解决方案:

#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

无论使用哪种方法，当我在教授的测试函数中为 log10() 输入 n = 10000 时，我得到这个错误:(异常(exception):当我的牛顿方法函数被使用时，log10() 偏离了，而这种二进制搜索方法在达到输入阈值之前相对准确，无论哪种方式，当 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')

这是测试函数:

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 ****'))

测试是这样设置的:

#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()

发布的其他解决方案似乎在运行全套测试输入时没有问题(我尽量不看发布的解决方案)。对此错误有任何见解吗？

---

似乎共识是数字的大小，但是，我教授的代码似乎适用于所有情况:

#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() 

结果:

     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

这其实是你理解数学而不是程序的问题。该算法很好，但提供的初始条件不是。

你这样定义inverse(f, delta):

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

所以您猜测 1000 = 10^x 的结果是 500.0，但是 10⁵⁰⁰ 肯定太大了。初始猜测应选择对 f 有效，而不是对 f 的逆选择。

我建议您使用 1 的猜测值进行初始化，即将该行替换为

guess = 1

它应该可以正常工作。

---

顺便说一句，你的二分搜索的初始条件也是错误的，因为你假设解决方案在 0 和 y 之间:

low, high = 0, float(y)

这对你的测试用例来说是正确的，但很容易构造反例，例如log₁₀ 0.1 (= -1), √0.36 (= 0.6) 等(你教授的find_bounds方法确实解决了√0.36的问题，但还是处理不了log₁₀ 0.1 问题。)