Search code examples
pythonrecursionmathformularecurrence

Refactor left recursion infinite loop

When you trying to implement this formula

enter image description here

in a straightforward manner using this code

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)

you will definitely end up with this error.

RecursionError: maximum recursion depth exceeded

So how to properly implement such class of problems using code? Is there a technique used to refactor this and get a solution? Or the only thing is to solve that mathematically first and only then feed into computation?

Solution

  • To expand on @Prune's suggestions, this recurrence can be written as two other recurrences, by moving either of the terms on the RHS to LHS:

    enter image description here

    These are non-cyclic because each term only depends on terms either before or after it, not both.

    However, they cannot be solved independently, because that requires knowledge of P1-l and Pl-1, where l = limit.

    Making some substitutions for convenience:

    enter image description here

    Note that, since the relations only contain linear terms (no e.g. Pn × Pn-1 etc.), all terms in their expansions will also be linear in the variables x, y:

    enter image description here

    These newly introduced coefficients also obey similar recurrence relations:

    enter image description here

    Their boundary conditions are:

    enter image description here

    Now that the boundary conditions are known, the coefficients can be computed iteratively. It only remains to equate the expressions for P1-l and Pl-1 from both recurrences:

    enter image description here

    Solving this pair of simultaneous equations gives P1-l and Pl-1, from which any Pn can be calculated (using the coefficients computed earlier):

    enter image description here

    Update: sample Python implementation

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

    Experimental results for 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   --

    Notes:

    • It is clear from the code and equations that the recurrence limits / boundary conditions can be arbitrary.
    • Numerical precision can be improved with summation techniques such as Kahan-Neumaier.
    • There also exists an alternative explicit formula for the iterative method used in solve_recurrence, which may produce more accurate results if used in conjunction with library functions.