https://www.v2ex.com/t/697506

main.py

import numpy as np
from numpy import matmul
from numpy.linalg import inv


class V2_697506(dict):
    def __missing__(self, key):
        x, y = key
        if y == 0:  # rule 1
            self[x, y] = 0
            return 0
        elif x == 0:  # rule 2
            self[x, y] = y
            return y
        elif x >= y:  # rule 3
            self[x, y] = self[y - 1, y]
            return self[x, y]
        else:  # rule 4
            self.solve_g(x + y)
            return self[x, y]

    def solve_g(self, g):
        alpha = (g - 1) // 2
        p0 = self[0, g]
        pa1 = self[alpha + 1, g - alpha - 1]
        if alpha == 1:
            self[1, g - 1] = (self[2, g - 2] + (g - 1) * self[0, g]) / g
            return
        lhs_k = np.eye(alpha)
        for i in range(alpha - 1):
            lhs_k[i][i + 1] = -(i + 1) / g
            lhs_k[i + 1][i] = -(g - i - 2) / g
        rhs = np.array([(g - 1) / g * p0, *([0] * (alpha - 2)), alpha / g * pa1])
        ilhs_k = inv(lhs_k)
        p = matmul(inv(lhs_k), rhs)
        for i in range(alpha):
            self[i + 1, g - i - 1] = p[i]


def test():
    f = V2_697506()
    for y in range(100):
        for x in range(100):
            if y == 0:
                if f[x, y] != 0:
                    raise ValueError
            elif x == 0:
                if f[x, y] != y:
                    raise ValueError
            elif x >= y:
                if f[x, y] != f[y - 1, y]:
                    raise ValueError
            else:
                lhs = f[x, y]
                rhs = (x * (f[x + 1, y - 1]) + y * (f[x - 1, y + 1])) / (x + y)
                if abs(lhs - rhs) >= 0.000001:
                    raise ValueError


if __name__ == '__main__':
    test()