Buggy python submission (verdict: Time Limit Exceeded):
N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0,0,0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1
n=x+y*2+z*3


dp=[[[-1 for _ in range(n+1)] for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0]=0
def calc(h,i,j):
    global n,N
    if h+i+j==0 or h<0 or h>n or i<0 or i>n or j<0:
        return 0
    if dp[h][i][j]!=-1:
        return dp[h][i][j]
    dp[h][i][j]=(calc(h-1,i,j)*h/N+calc(h+1,i-1,j)*i/N+calc(h,i+1,j-1)*j/N+1)*N/(h+i+j)
    return dp[h][i][j]

a=calc(x,y,z)
print(dp[x][y][z])


Fixed submission:
import numpy as np
from numba import njit

N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0,0,0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1


dp=np.full((N+1,N+1,N+1),-1,dtype=np.float64)

dp[0][0][0]=0

@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def calc(h,i,j,N,dp):
    if h+i+j==0:
        return 0
    if dp[h][i][j]>=0:
        return dp[h][i][j]
    A=1
    if h>=1:
        A+=calc(h-1,i,j,N,dp)*h/N
    if i>=1:
        A+=calc(h+1,i-1,j,N,dp)*i/N
    if j>=1:
        A+=calc(h,i+1,j-1,N,dp)*j/N
    A=A*N/(h+i+j)
    dp[h][i][j]=A
    return A

a=calc(x,y,z,N,dp)

print(dp[x][y][z])

Buggy python submission (verdict: Runtime Error):
N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1
n=x+y*2+z*3


dp=[[[-1 for _ in range(n+1)] for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0]=0
def calc(h,i,j):
    global n,N
    if h+i+j==0 or h<0 or h>n or i<0 or i>n or j<0:
        return 0
    if dp[h][i][j]!=-1:
        return dp[h][i][j]
    dp[h][i][j]=(calc(h-1,i,j)*h/N+calc(h+1,i-1,j)*i/N+calc(h,i+1,j-1)*j/N+1)*N/(h+i+j)
    return dp[h][i][j]

a=calc(x,y,z)
print(dp[x][y][z])

Fixed submission:
import numpy as np
from numba import njit

N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0,0,0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1


dp=np.full((N+1,N+1,N+1),-1,dtype=np.float64)

dp[0][0][0]=0

@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def calc(h,i,j,N,dp):
    if h+i+j==0:
        return 0
    if dp[h][i][j]>=0:
        return dp[h][i][j]
    A=1
    if h>=1:
        A+=calc(h-1,i,j,N,dp)*h/N
    if i>=1:
        A+=calc(h+1,i-1,j,N,dp)*i/N
    if j>=1:
        A+=calc(h,i+1,j-1,N,dp)*j/N
    A=A*N/(h+i+j)
    dp[h][i][j]=A
    return A

a=calc(x,y,z,N,dp)

print(dp[x][y][z])

Buggy python submission (verdict: Time Limit Exceeded):
N=int(input())
a=[int(i) for i in input().split()]

x,y,z=a.count(1),a.count(2),a.count(3)
n=x+y*2+z*3


dp=[[[-1 for _ in range(n+1)] for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0]=0
def calc(h,i,j):
    global n,N
    if h+i+j==0 or h<0 or h>n or i<0 or i>n or j<0:
        return 0
    if dp[h][i][j]!=-1:
        return dp[h][i][j]
    dp[h][i][j]=(calc(h-1,i,j)*h/N+calc(h+1,i-1,j)*i/N+calc(h,i+1,j-1)*j/N+1)*N/(h+i+j)
    return dp[h][i][j]

a=calc(x,y,z)
print(dp[x][y][z])


Fixed submission:
import numpy as np
from numba import njit

N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0,0,0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1


dp=np.full((N+1,N+1,N+1),-1,dtype=np.float64)

dp[0][0][0]=0

@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def calc(h,i,j,N,dp):
    if h+i+j==0:
        return 0
    if dp[h][i][j]>=0:
        return dp[h][i][j]
    A=1
    if h>=1:
        A+=calc(h-1,i,j,N,dp)*h/N
    if i>=1:
        A+=calc(h+1,i-1,j,N,dp)*i/N
    if j>=1:
        A+=calc(h,i+1,j-1,N,dp)*j/N
    A=A*N/(h+i+j)
    dp[h][i][j]=A
    return A

a=calc(x,y,z,N,dp)

print(dp[x][y][z])

Buggy python submission (verdict: Time Limit Exceeded):
N=int(input())
a=[int(i) for i in input().split()]

x,y,z=a.count(1),a.count(2),a.count(3)
n=x+y*2+z*3


dp=[[[0 for _ in range(n+2)] for _ in range(n+2)] for _ in range(n+2)]

for j in range(n+1):
    for i in range(n+1):
        for h in range(n+1):
            if h+i+j==0 or h+i*2+j*3>n:
                continue
            dp[h][i][j]=(dp[h-1][i][j]*h/N+dp[h+1][i-1][j]*i/N+dp[h][i+1][j-1]*j/N+1)*N/(h+i+j)
            if h==x and i==y and j==z:
                print(dp[x][y][z])
                exit()


Fixed submission:
import numpy as np
from numba import njit

N=int(input())
a=[int(i) for i in input().split()]
x,y,z=0,0,0
for i in a:
    if i==1:
        x+=1
    elif i==2:
        y+=1
    elif i==3:
        z+=1


dp=np.full((N+1,N+1,N+1),-1,dtype=np.float64)

dp[0][0][0]=0

@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def calc(h,i,j,N,dp):
    if h+i+j==0:
        return 0
    if dp[h][i][j]>=0:
        return dp[h][i][j]
    A=1
    if h>=1:
        A+=calc(h-1,i,j,N,dp)*h/N
    if i>=1:
        A+=calc(h+1,i-1,j,N,dp)*i/N
    if j>=1:
        A+=calc(h,i+1,j-1,N,dp)*j/N
    A=A*N/(h+i+j)
    dp[h][i][j]=A
    return A

a=calc(x,y,z,N,dp)

print(dp[x][y][z])