重构左递归死循环

Question

当你尝试实现这个公式时

直接使用此代码

def P(n, p, limit):
    if n == limit:
        return 1

    if n == -limit:
        return 0

    return p*P(n+1, p, limit) + (1-p)*P(n-1, p, limit)

你肯定会遇到这个错误。

RecursionError: maximum recursion depth exceeded

那么如何正确地使用代码实现这样的class问题呢？有没有一种技术可以用来重构它并获得解决方案？或者唯一的办法就是先用数学方法解决这个问题，然后再输入计算？

Answer 1

递归取决于将问题分解成多个部分，每个部分都更接近于一个基本情况。您给出的公式是无限的，因为它取决于同样未定义的级数项。看似简单的数学符号并不能保证可计算性。

在这种情况下，您必须以非循环方式重写递归关系。你仍然可以使用递归，但必须有一个明确的结束。在给定的公式中，P[n]依赖于P[n+1]和P[n-1]，每个直接依赖于P[n]。这是算法循环。

您需要从基本案例中提取算法，从末端开始，逐步深入。此递归关系的不幸结果之一是您必须求解整个范围 [-limit, limit]为了得到你想要的职位的价值。基本案例值从两端波动到您的观点。

没有工具可以象征性地为您解开这个谜团。最重要的是，请注意您的基本情况不会出现在等式中。即使您将这些作为附加方程式提供，从这种关系中导出指数级数也超出了当今求解器的分析能力。

Answer 2

为了扩展 @Prune 的建议，通过将 RHS 上的任一术语移动到 LHS，可以将此循环写成两个其他循环：

这些是非循环的，因为每个项仅取决于或之前的项，而不是两者。

然而，它们不能独立求解，因为这需要 P_1-l 和 P_l-1 的知识，其中 l = limit.

为方便起见做一些替换：

请注意，由于关系仅包含线性项（没有例如 P_n × P_n-1 等），其展开式中的所有项在变量 x, y:

中也将是线性的

这些新引入的系数也服从类似的递归关系：

它们的边界条件是：

现在已知边界条件，可以迭代计算系数。剩下的只是将 P_1-l 和 P_l-1 的表达式从两个递归中等同起来：

求解这对联立方程得到P_1-l和P_l-1，从中任意P_n可以计算（使用之前计算的系数）：

更新：样本Python实施

# solve a recurrence of the form
# Tn = c1 * Tn-1 + c2 * Tn-2, n >= 0
def solve_recurrence(c1, c2, t1, t2, n):
    if n <= 0: return t2
    if n == 1: return t1
    a, b = t2, t1
    for i in range(1, n):
        a, b = b, c1 * b + c2 * a
    return b

# solve the loop recurrence
def solve_loop(limit, p, n):
    assert(limit > 0 and abs(n) <= limit)
    assert(p > 0 and p < 1)

    # corner cases
    if n == -limit: return 0
    if n ==  limit: return 1

    # constants
    a, c = 1 / p, 1 / (1 - p)
    b, d = 1 - a, 1 - c

    # coefficients at 2l - 1
    e = 2 * limit - 1
    l = solve_recurrence(a, b, 1, 0, e)
    m = solve_recurrence(a, b, 0, 0, e)
    s = solve_recurrence(c, d, 1, 0, e)
    t = solve_recurrence(c, d, 0, 1, e)
    # note that m can just be 0 but was left in for generality

    # solving the simultaneous equations
    # to get P at 1 - l and l - 1
    x = (s * m + t) / (1 - s * l)
    y = (l * t + m) / (1 - l * s)

    # iterate from the closest side to improve precision
    if n < 0:
        return solve_recurrence(a, b, x, 0, limit + n)
    else:
        return solve_recurrence(c, d, y, 1, limit - n)

limit = 10, p = 0.1 的实验结果：

n       | P(n)            abs. error
========================================
-10     | 0.0000000E+00   --
-9      | 6.5802107E-19   0.0000000E+00
-8      | 6.5802107E-18   1.7995889E-30
-7      | 5.9879917E-17   1.7995889E-29
-6      | 5.3957728E-16   5.0003921E-28
-5      | 4.8568535E-15   8.1994202E-27
-4      | 4.3712339E-14   8.9993252E-27
-3      | 3.9341171E-13   5.1002066E-25
-2      | 3.5407061E-12   5.9938283E-25
-1      | 3.1866355E-11   5.9339980E-18
 0      | 2.8679726E-10   5.3406100E-17
 1      | 2.5811749E-09   3.3000371E-21
 2      | 2.3230573E-08   2.6999175E-20
 3      | 2.0907516E-07   2.2999591E-19
 4      | 1.8816764E-06   7.9981086E-19
 5      | 1.6935088E-05   1.9989978E-18
 6      | 1.5241579E-04   3.5000757E-16
 7      | 1.3717421E-03   1.5998487E-15
 8      | 1.2345679E-02   3.2000444E-14
 9      | 1.1111111E-01   6.9999562E-14
 10     | 1.0000000E+00   --

备注：

从代码和方程式可以清楚地看出，递归限制/边界条件可以是任意的。
可以使用 Kahan-Neumaier 等求和技术提高数值精度。
对于 solve_recurrence 中使用的迭代方法，还有一个替代显式公式，如果与库函数结合使用，可能会产生更准确的结果。

Answer 3

我有另一个使用符号计算的解决方案。

class X(object):
    """Representing unknown value which will be resolved in future"""
    def __init__(self, value=1.0):
        self.value = value

    def __add__(self, other):
        cls = self.__class__
        return cls(value=self.value+other.value)

    def __sub__(self, other):
        cls = self.__class__
        return cls(value=self.value-other.value)

    def __mul__(self, other):
        cls = self.__class__
        if isinstance(other, cls):
            return cls(value=self.value*other.value)
        elif isinstance(other, numbers.Real):
            return cls(value=self.value*other)
        else:
            raise TypeError(f'Cannot multiply with type {type(other)}')

    def __rmul__(self, other):
        return self.__mul__(other)

我们没有立即获得数字结果，而是依赖于这样一个事实，即结果值将以某种方式表示为

每个中间项将表示为

这就是我们为对象定义 add/subtract/multiplicate 操作的原因。

所以在代码中有这样一个符号的表示，我们可以用类似于公认答案的方法来求解这个方程（我们从 -limit 和 limit 两边出发，并将其他项表示为未知产品）。当到达基点（0 状态）时，我们将得到概率与该比率的数值之比，这使我们可以轻松计算结果

_memory = []
def fill_memory(n):
    global _memory
    _memory = [(0, X(1))]*(2*n+1)
    _memory[-n] = (0, X(0))
    _memory[n] = (1, X(0))


def resolve(n, p):
    q = 1-p
    for i in reversed(range(n)):
        # resolve from right to center
        prob_next = _memory[i+1]
        prob_prev = _memory[i-i]
        _memory[i] = (
            q*prob_prev[0] + p*prob_next[0],
            q*prob_prev[1] + p*prob_next[1],
        )

        # resolve from left to center
        prob_prev = _memory[-i-1]
        prob_next = _memory[-i+1]
        _memory[-i] = (
            q*prob_prev[0] + p*prob_next[0],
            q*prob_prev[1] + p*prob_next[1],
        )

    # calculate base probability
    solution = _memory[0]
    coef = (X(value=1) - solution[1]).value
    p0 = solution[0] / coef
    return p0

我用这种方法获得了正确且非常精确的结果（我使用的是 p=0.6 和 N=4（限制））

def main():
    N = 4
    p = 0.6
    fill_memory(N)
    p0 = resolve(N, p)
    print(f'Probability for N={N} steps is {p0*100:.4}%')
    for i in range(2*N+1):
        print(i-N, '->', _memory[i-N])

4 - 83.505%
3 - 77.143%
2 - 69.231%
1 - 60.0%

如您所见，这与数学求解的结果相同。

重构左递归死循环

Refactor left recursion infinite loop

python

recursion

math

formula

recurrence