Buggy python submission (verdict: Time Limit Exceeded):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

N = int(readline())
A = np.array(read().split(), np.int8)

counter = np.bincount(A,minlength=4)
S = A.sum()

B = (counter[2]+counter[3])
C = counter[3]

dp = np.zeros((B+1,C+1),np.float64)

p2 = np.tile(np.arange(B+1)[:,None],(1,C+1))/N
p3 = np.tile(np.arange(C+1)[None,:],(B+1,1))/N
p_rest = 1-p2-p3

bc = np.arange(B+1)[:,None]*2+np.arange(C+1)[None,:]*3

for n in range(1,S+1):
    prev = dp
    one = n-bc
    p1 = np.maximum(0,one/N)
    dp = np.ones_like(prev)
    dp += p1*prev
    dp[1:] += p2[1:]*prev[:-1]
    dp[:-1,1:] += p3[:-1,1:]*prev[1:,:-1]
    dp /= (p1+p2+p3)

answer = dp[B-C,C]
print(answer)

Fixed submission:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
N = int(input())

sushi = np.bincount(np.array(input().split(), dtype=np.int8),minlength=4)
total = (sushi * np.arange(4)).sum()

# a,b,c = 1個、2個、3個
# S = a + 2b + 3c：total スシ
# S の小さい方から求めていく
# S ごとには、b,c だけ持てば a が分かる

b = np.arange(N+1)[:,None]
c = np.arange(N+1)[None,:]
S_bc = 2 * b + 3 * c
dp = np.zeros((N+1,N+1))
for S in range(1,total+1):
    prev = dp
    a = S - S_bc
    dp = prev * a + N
    dp[1:] += prev[:-1] * b[1:]
    dp[:-1,1:] += prev[1:,:-1] * c[:,1:]
    dp /= (a+b+c)
    dp[a < 0] = 0

answer = dp[sushi[2],sushi[3]]
print(answer)

Buggy python submission (verdict: Runtime Error):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

N,*A = map(int,read().split())

from functools import lru_cache
from collections import Counter

@lru_cache(None)
def F(one,two,thr):
    if one==two==thr==0:
        return 0
    x = 1
    if one:
        x += one/N*F(one-1,two,thr)
    if two:
        x += two/N*F(one+1,two-1,thr)
    if thr:
        x += thr/N*F(one,two+1,thr-1)
    x *= N; x /= (one+two+thr)
    return x

counter = Counter(A)

answer = F(*[counter[n] for n in [1,2,3]])
print(answer)

Fixed submission:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
N = int(input())

sushi = np.bincount(np.array(input().split(), dtype=np.int8),minlength=4)
total = (sushi * np.arange(4)).sum()

# a,b,c = 1個、2個、3個
# S = a + 2b + 3c：total スシ
# S の小さい方から求めていく
# S ごとには、b,c だけ持てば a が分かる

b = np.arange(N+1)[:,None]
c = np.arange(N+1)[None,:]
S_bc = 2 * b + 3 * c
dp = np.zeros((N+1,N+1))
for S in range(1,total+1):
    prev = dp
    a = S - S_bc
    dp = prev * a + N
    dp[1:] += prev[:-1] * b[1:]
    dp[:-1,1:] += prev[1:,:-1] * c[:,1:]
    dp /= (a+b+c)
    dp[a < 0] = 0

answer = dp[sushi[2],sushi[3]]
print(answer)

Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
N = int(input())

sushi = np.bincount(np.array(input().split(), dtype=np.int8),minlength=4)
total = (sushi * np.arange(4)).sum()

# a,b,c = 1個、2個、3個
# S = a + 2b + 3c：total スシ
# S の小さい方から求めていく
# S ごとには、b,c だけ持てば a が分かる

b = np.arange(N+1)[:,None]
c = np.arange(N+1)[None,:]
S_bc = 2 * b + 3 * c
dp = np.zeros((N+1,N+1))
for S in range(1,total+1):
    prev = dp
    a = S - S_bc
    dp = np.ones_like(prev) * N
    dp += prev * a# 1個の皿
    dp[1:] += prev[:-1] * b[1:] # 2個の皿
    dp[:-1,1:] += prev[1:,:-1] * c[:,1:] # 3個の皿
    dp /= (a+b+c)
    dp[a < 0] = 0

answer = dp[sushi[2],sushi[3]]
print(answer)

Fixed submission:
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

import numpy as np
N = int(input())

sushi = np.bincount(np.array(input().split(), dtype=np.int8),minlength=4)
total = (sushi * np.arange(4)).sum()

# a,b,c = 1個、2個、3個
# S = a + 2b + 3c：total スシ
# S の小さい方から求めていく
# S ごとには、b,c だけ持てば a が分かる

b = np.arange(N+1)[:,None]
c = np.arange(N+1)[None,:]
S_bc = 2 * b + 3 * c
dp = np.zeros((N+1,N+1))
for S in range(1,total+1):
    prev = dp
    a = S - S_bc
    dp = prev * a + N
    dp[1:] += prev[:-1] * b[1:]
    dp[:-1,1:] += prev[1:,:-1] * c[:,1:]
    dp /= (a+b+c)
    dp[a < 0] = 0

answer = dp[sushi[2],sushi[3]]
print(answer)