Buggy python submission (verdict: Time Limit Exceeded):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]

    z = r[0]
    for j in range(E):
        x = z[j]
        if j > 0:
            m = z[j-1]
        for i in range(F-j):
            if not(i or j):
                continue
            w = n
            if i > 0:
                w += x[i-1] * i
            if j > 0:
                w += m[i+1] * j
            x[i] = w / (i+j)

    for k in range(1, D):
        y,z = z,r[k]
        G = F - k

        x = z[0]
        p = y[1]
        for i in range(G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j,p,t in zip(range(1,E-k), y[2:],z[1:]):
            for i,w,v in zip(range(G-j), p, x[1:]):
                t[i] = (n + t[i-1] * i + v * j + w * k) / (i+j+k)
            x = t

    return "{:0.9f}".format(r[C][B][A])


print(main())


Fixed submission:
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())

Buggy python submission (verdict: Time Limit Exceeded):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())


Fixed submission:
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())

Buggy python submission (verdict: Time Limit Exceeded):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(D):
        t0 = r[k]
        G = F - k
        for j in range(E-k):
            tj = t0[j]
            if j > 0:
                tm = t0[j-1]
            if k > 0:
                tp = t1[j+1]
            for i in range(G-j):
                if not(i or j or k):
                    continue
                tr = n
                if i > 0:
                    tr += tj[i-1] * i
                if j > 0:
                    tr += tm[i+1] * j
                if k > 0:
                    tr += tp[i] * k
                tj[i] = tr / (i+j+k)
        t1 = t0


    return "{:0.9f}".format(r[C][B][A])


print(main())


Fixed submission:
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())

Buggy python submission (verdict: Time Limit Exceeded):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 4
    for c in a:
        t[c] += 1
    s = [0] * 4
    s[3] = t[3]
    for i in range(2,-1,-1):
        s[i] = s[i+1] + t[i]

    r = [[[0]*(n+1) for _ in range(n+1)] for __ in range(n+1)]
    for c in range(1,n+1):
        for i in range(c,-1,-1):
            if c - i > s[2]:
                break
            ri = r[i]
            if i > 0:
                rm = r[i-1]
            if i < n:
                rp = r[i+1]
            for j in range(c-i,-1,-1):
                k = c - i - j
                if k > s[3]:
                    break
                tr = n / c
                if i > 0:
                    tr += rm[j][k] * i / c
                if j > 0:
                    tr += rp[j-1][k] * j / c
                if k > 0:
                    tr += ri[j+1][k-1] * k / c
                ri[j][k] = tr


    return "{:0.9f}".format(r[t[1]][t[2]][t[3]])


print(main())


Fixed submission:
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())

Buggy python submission (verdict: Time Limit Exceeded):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 4
    for c in a:
        t[c] += 1

    fm = {}
    fm[(n,0,0,0)] = 0
    def f(t):
        if t in fm:
            return fm[t]
        c = n-t[0]
        r = 1 / (c / n)
        ta = list(t)
        for i in range(1,4):
            if t[i] == 0:
                continue
            ta[i-1] += 1
            ta[i] -= 1
            u = f(tuple(ta))
            r += u * t[i] / c
            ta[i-1] -= 1
            ta[i] += 1
        fm[t] = r
        return r
    r = f(tuple(t))

    return "{:0.9f}".format(r)


print(main())


Fixed submission:
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = LI()
    t = [0] * 3
    for c in a:
        t[c-1] += 1
    A,B,C = t
    D = C + 1
    E = B + D
    F = A + E

    r = [[[0]*(F-k-j) for j in range(E-k)] for k in range(D)]
    for k in range(1):
        z = r[k]
        G = F - k
        for j in range(E-k):
            x = z[j]
            if j > 0:
                m = z[j-1]
            l = k + j
            for i in range(G-j):
                if not(i or l):
                    continue
                w = n
                if i > 0:
                    w += x[i-1] * i
                if j > 0:
                    w += m[i+1] * j
                x[i] = w / (i+l)


    y = r[0]
    for k in range(1, D):
        z = r[k]
        G = F - k

        x = z[0]
        p = y[1]
        x[0] = (n + p[0] * k) / k

        for i in range(1,G):
            x[i] = (n + x[i-1] * i + p[i] * k) / (i+k)

        for j in range(1,E-k):
            x = z[j]
            m = z[j-1]
            p = y[j+1]
            l = j + k
            x[0] = (n + m[1] * j + p[0] * k) / l

            for i in range(1,G-j):
                x[i] = (n + x[i-1] * i + m[i+1] * j + p[i] * k) / (i+l)

        y = z


    return "{:0.9f}".format(r[C][B][A])


print(main())