Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    for i in range(n + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(n + 1):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(n + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                d[i][j][k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    d = defaultdict(int)
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    for i in range(n+1):
        for j in range(n+1):
            for k in range(n+1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * d[(k, j + 1, i - 1)]
                if j: res += j * d[(k + 1, j - 1, i)]
                if k: res += k * d[(k - 1, j, i)]
                d[(k, j, i)] = res / l
    print(d[(c1, c2, c3)])
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    d = defaultdict(int)
    @lru_cache(maxsize=300 ** 3)
    def f(c1, c2, c3):
        tmp = 1 / (1 - (n - c1 - c2 - c3) / n) if c1 != 0 or c2 != 0 or c3 != 0 else 0
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c2 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c3 / n / (1 - (n - c1 - c2 - c3) / n)
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    d = defaultdict(int)
    @lru_cache(maxsize=300 ** 3)
    def f(c1, c2, c3):
        tmp = 1 / (1 - (n - c1 - c2 - c3) / n) if c1 != 0 or c2 != 0 or c3 != 0 else 0
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c2 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c3 / n / (1 - (n - c1 - c2 - c3) / n)
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10

#solve
def solve():
    n = II()
    d = defaultdict(int)
    def f(c1, c2, c3):
        if d[(c1,c2,c3)]: return d[(c1,c2,c3)]
        tmp = 1 / (1 - (n - c1 - c2 - c3) / n) if c1 != 0 or c2 != 0 or c3 != 0 else 0
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c2 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c3 / n / (1 - (n - c1 - c2 - c3) / n)
        d[(c1, c2, c3)] = tmp
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Time Limit Exceeded):
#!/usr/bin/env python3
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10

#solve
def solve():
    n = II()
    d = [[[0] * (n+1) for _ in range(n+1)] for _ in range(n+1)]
    def f(c1, c2, c3):
        if d[c1][c2][c3]: return d[c1][c2][c3]
        tmp = 1 / (1 - (n - c1 - c2 - c3) / n) if c1 != 0 or c2 != 0 or c3 != 0 else 0
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c2 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c3 / n / (1 - (n - c1 - c2 - c3) / n)
        d[c1][c2][c3] = tmp
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Runtime Error):
#!/usr/bin/env python3
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10

#solve
def solve():
    n = II()
    d = [[[0] * n for _ in range(n)] for _ in range(n)]
    def f(c1, c2, c3):
        if d[c1][c2][c3]: return d[c1][c2][c3]
        tmp = 1 / (1 - (n - c1 - c2 - c3)/n)
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c2 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c3 / n / (1 - (n - c1 - c2 - c3) / n)
        d[c1][c2][c3] = tmp
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()

Buggy python submission (verdict: Runtime Error):
#!/usr/bin/env python3
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10

#solve
def solve():
    n = II()
    d = [[[0] * n for _ in range(n)] for _ in range(n)]
    def f(c1, c2, c3):
        if d[c1][c2][c3]: return d[c1][c2][c3]
        tmp = 1 / (1 - (n - c1 - c2 - c3)/n)
        if c1: tmp += f(c1 - 1, c2, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c2: tmp += f(c1 + 1, c2 - 1, c3) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        if c3: tmp += f(c1, c2 + 1, c3 - 1) * c1 / n / (1 - (n - c1 - c2 - c3) / n)
        d[c1][c2][c3] = tmp
        return tmp
    c = LI()
    print(f(c.count(1), c.count(2), c.count(3)))
    return


#main
if __name__ == '__main__':
    solve()


Fixed submission:
#!/usr/bin/env python3
from functools import lru_cache
from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def LS(): return list(map(list, input().split()))
def S(): return list(input().rstrip())
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = 1e10
#solve
def solve():
    n = II()
    c = LI()
    c1, c2, c3 = [c.count(i) for i in range(1, 4)]
    d = [[[0] * (n + 2) for _ in range(n + 2)] for _ in range(n + 2)]
    for i in range(c3 + 1):
        di = d[i]
        dim1 = d[i - 1]
        for j in range(c3 + c2 + 1 - i):
            dim1jp1 = dim1[j + 1]
            dijm1 = di[j - 1]
            dij = di[j]
            for k in range(c3 + c2 + c1 - i - j + 1):
                l = i + j + k
                if l == 0: continue
                res = n
                if i: res += i * dim1jp1[k]
                if j: res += j * dijm1[k + 1]
                if k: res += k * dij[k-1]
                dij[k] = res / l
    print(d[c3][c2][c1])
    return


#main
if __name__ == '__main__':
    solve()