Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
from collections import Counter, deque
from collections import defaultdict
from itertools import combinations, permutations, accumulate, groupby, product
from bisect import bisect_left,bisect_right
from heapq import heapify, heappop, heappush
from math import floor, ceil,pi,factorial
from operator import itemgetter
def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))
def LI2(): return [int(input()) for i in range(n)]
def MXI(): return [[LI()]for i in range(n)]
def SI(): return input().rstrip()
def printns(x): print('\n'.join(x))
def printni(x): print('\n'.join(list(map(str,x))))
inf = 10**17
mod = 10**9 + 7

n=I()
lis=LI()
c=Counter(lis)
dp=[[[-1]*(n+1) for i in range(n+1)] for j in range(n+1)]
dp[0][0][0]=0
def dfs(x,y,z):
    if x+y+z>n or (x<0 or y<0 or z<0) or (x>n or y>n or z>n):
        return 0
    elif dp[x][y][z]>=0:
        return dp[x][y][z]
    else:
        dp[x][y][z]=(n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
        return (n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
print(dfs(c[1],c[2],c[3]))


Fixed submission:
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    # 0個の皿がx0枚、1個の皿がx1枚、2個の皿がx2枚、3個の皿がx3枚
    # x0=n-x1-x2-x3
    ca={1:a.count(1),2:a.count(2),3:a.count(3)}
    memo=[[[-1]*(n+2) for _ in range(n+2)] for _ in range(n+2)]
    memo[0][0][0]=0
    for x3 in range(ca[3]+1):
        for x2 in range(ca[2]+ca[3]+1-x3):
            for x1 in range(ca[1]+ca[2]+ca[3]+1-x3-x2):
                if (x1,x2,x3)==(0,0,0):continue
                tmp=0
                if x1>0:tmp+=memo[x1-1][x2][x3]*x1
                if x2>0:tmp+=memo[x1+1][x2-1][x3]*x2
                if x3>0:tmp+=memo[x1][x2+1][x3-1]*x3
                tmp+=n
                tmp/=x1+x2+x3
                memo[x1][x2][x3]=tmp
    print(memo[ca[1]][ca[2]][ca[3]])
if __name__=='__main__':
    main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
from collections import Counter, deque
from collections import defaultdict
from itertools import combinations, permutations, accumulate, groupby, product
from bisect import bisect_left,bisect_right
from heapq import heapify, heappop, heappush
from math import floor, ceil,pi,factorial
from operator import itemgetter
def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))
def LI2(): return [int(input()) for i in range(n)]
def MXI(): return [[LI()]for i in range(n)]
def SI(): return input().rstrip()
def printns(x): print('\n'.join(x))
def printni(x): print('\n'.join(list(map(str,x))))
inf = 10**17
mod = 10**9 + 7

n=I()
lis=LI()
c=Counter(lis)
dp=[[[-1]*(n+1) for i in range(n+1)] for j in range(n+1)]
dp[0][0][0]=0
def dfs(x,y,z):
    if x+y+z>n or (x<0 or y<0 or z<0) or (x>n or y>n or z>n):
        return 0
    elif dp[x][y][z]>=0:
        return dp[x][y][z]
    else:
        dp[x][y][z]=(n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
        return (n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
print(dfs(c[1],c[2],c[3]))


Fixed submission:
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    # 0個の皿がx0枚、1個の皿がx1枚、2個の皿がx2枚、3個の皿がx3枚
    # x0=n-x1-x2-x3
    ca={1:a.count(1),2:a.count(2),3:a.count(3)}
    memo=[[[-1]*(n+2) for _ in range(n+2)] for _ in range(n+2)]
    memo[0][0][0]=0
    for x3 in range(ca[3]+1):
        for x2 in range(ca[2]+ca[3]+1-x3):
            for x1 in range(ca[1]+ca[2]+ca[3]+1-x3-x2):
                if (x1,x2,x3)==(0,0,0):continue
                tmp=0
                if x1>0:tmp+=memo[x1-1][x2][x3]*x1
                if x2>0:tmp+=memo[x1+1][x2-1][x3]*x2
                if x3>0:tmp+=memo[x1][x2+1][x3-1]*x3
                tmp+=n
                tmp/=x1+x2+x3
                memo[x1][x2][x3]=tmp
    print(memo[ca[1]][ca[2]][ca[3]])
if __name__=='__main__':
    main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
from collections import Counter, deque
from collections import defaultdict
from itertools import combinations, permutations, accumulate, groupby, product
from bisect import bisect_left,bisect_right
from heapq import heapify, heappop, heappush
from math import floor, ceil,pi,factorial
from operator import itemgetter
def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))
def LI2(): return [int(input()) for i in range(n)]
def MXI(): return [[LI()]for i in range(n)]
def SI(): return input().rstrip()
def printns(x): print('\n'.join(x))
def printni(x): print('\n'.join(list(map(str,x))))
inf = 10**17
mod = 10**9 + 7

n=I()
lis=LI()
c=Counter(lis)
dp=[[[-1]*(n+1) for i in range(n+1)] for j in range(n+1)]
dp[0][0][0]=0
def dfs(x,y,z):
    if x+y+z>n or (x<0 or y<0 or z<0) or (x>n or y>n or z>n):
        return 0
    elif dp[x][y][z]>=0:
        return dp[x][y][z]
    else:
        dp[x][y][z]=(n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
        return (n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
print(dfs(c[1],c[2],c[3]))
'''print(dp[0])
print(dp[1])
print(dp[2])
print(dp[3])'''
#print(dp[c[1]][c[2]][c[3]])

Fixed submission:
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    # 0個の皿がx0枚、1個の皿がx1枚、2個の皿がx2枚、3個の皿がx3枚
    # x0=n-x1-x2-x3
    ca={1:a.count(1),2:a.count(2),3:a.count(3)}
    memo=[[[-1]*(n+2) for _ in range(n+2)] for _ in range(n+2)]
    memo[0][0][0]=0
    for x3 in range(ca[3]+1):
        for x2 in range(ca[2]+ca[3]+1-x3):
            for x1 in range(ca[1]+ca[2]+ca[3]+1-x3-x2):
                if (x1,x2,x3)==(0,0,0):continue
                tmp=0
                if x1>0:tmp+=memo[x1-1][x2][x3]*x1
                if x2>0:tmp+=memo[x1+1][x2-1][x3]*x2
                if x3>0:tmp+=memo[x1][x2+1][x3-1]*x3
                tmp+=n
                tmp/=x1+x2+x3
                memo[x1][x2][x3]=tmp
    print(memo[ca[1]][ca[2]][ca[3]])
if __name__=='__main__':
    main()

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
from collections import Counter, deque
from collections import defaultdict
from itertools import combinations, permutations, accumulate, groupby, product
from bisect import bisect_left,bisect_right
from heapq import heapify, heappop, heappush
from math import floor, ceil,pi,factorial
from operator import itemgetter
def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))
def LI2(): return [int(input()) for i in range(n)]
def MXI(): return [[LI()]for i in range(n)]
def SI(): return input().rstrip()
def printns(x): print('\n'.join(x))
def printni(x): print('\n'.join(list(map(str,x))))
inf = 10**17
mod = 10**9 + 7

n=I()
lis=LI()
c=Counter(lis)
dp=[[[-1]*(n+1) for i in range(n+1)] for j in range(n+1)]
dp[0][0][0]=0
def dfs(x,y,z):
    if x+y+z>n or (x<0 or y<0 or z<0) or (x>n or y>n or z>n):
        return 0
    elif dp[x][y][z]>=0:
        return dp[x][y][z]
    else:
        dp[x][y][z]=(n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
        return (n+x*dfs(x-1,y,z)+y*dfs(x+1,y-1,z)+z*dfs(x,y+1,z-1))/(x+y+z)
for i in range(n+1):
    for j in range(n+1):
        for k in range(n+1):
            dfs(i,j,k)
'''print(dp[0])
print(dp[1])
print(dp[2])
print(dp[3])'''
print(dp[c[1]][c[2]][c[3]])

Fixed submission:
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    # 0個の皿がx0枚、1個の皿がx1枚、2個の皿がx2枚、3個の皿がx3枚
    # x0=n-x1-x2-x3
    ca={1:a.count(1),2:a.count(2),3:a.count(3)}
    memo=[[[-1]*(n+2) for _ in range(n+2)] for _ in range(n+2)]
    memo[0][0][0]=0
    for x3 in range(ca[3]+1):
        for x2 in range(ca[2]+ca[3]+1-x3):
            for x1 in range(ca[1]+ca[2]+ca[3]+1-x3-x2):
                if (x1,x2,x3)==(0,0,0):continue
                tmp=0
                if x1>0:tmp+=memo[x1-1][x2][x3]*x1
                if x2>0:tmp+=memo[x1+1][x2-1][x3]*x2
                if x3>0:tmp+=memo[x1][x2+1][x3-1]*x3
                tmp+=n
                tmp/=x1+x2+x3
                memo[x1][x2][x3]=tmp
    print(memo[ca[1]][ca[2]][ca[3]])
if __name__=='__main__':
    main()