Buggy python submission (verdict: Time Limit Exceeded):
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])

Fixed submission:
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])

Buggy python submission (verdict: Time Limit Exceeded):
N = int(input())
A = input().split()

px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nz+1) for y in range(ny+1)] for x in range(nx+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[x-1][y][z]*x/s
        if y > 0:
          t += dp[x+1][y-1][z]*y/s
        if z > 0:
          t += dp[x][y+1][z-1]*z/s
        
        dp[x][y][z] += t


print(dp[px][py][pz])

Fixed submission:
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])

Buggy python submission (verdict: Time Limit Exceeded):
N = int(input())
A = list(map(int, input().split()))
pos = [0, 0, 0]
for a in A:
  if a == 1:
    pos[0] += 1
  elif a == 2:
    pos[1] += 1
  elif a == 3:
    pos[2] += 1


dp = [[[0]*(N+1) for y in range(N+1)] for x in range(N+1)]
transition = [[1, 0, 0, 0], [-1, 1, 0, 1], [0, -1, 1, 2]]

def get_expectation(pos):
  global dp
  px, py, pz = pos
  if px == py == pz == 0:
    return 0
  elif px<0 or px>N or py<0 or py>N or pz<0 or pz>N:
    return 0
  else:
    if dp[px][py][pz] > 0:
      return dp[px][py][pz]
    else:
      for dx, dy, dz, idx in transition:
        s = px + py + pz
        dp[px][py][pz] += (get_expectation([px-dx, py-dy, pz-dz]) + N/s)*[px, py, pz][idx]/s
      
      return dp[px][py][pz]


print(get_expectation(pos))

Fixed submission:
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])

Buggy python submission (verdict: Time Limit Exceeded):
N = int(input())
A = list(map(int, input().split()))
pos = [0, 0, 0]
for a in A:
  if a == 1:
    pos[0] += 1
  elif a == 2:
    pos[1] += 1
  elif a == 3:
    pos[2] += 1


dp = [[[0]*(N+1) for y in range(N+1)] for x in range(N+1)]


def get_expectation(pos):
  global dp
  px, py, pz = pos
  if px == py == pz == 0:
    return 0
  elif px<0 or px>N or py<0 or py>N or pz<0 or pz>N:
    return 0
  else:
    if dp[px][py][pz] > 0:
      return dp[px][py][pz]
    else:
      transition = [[1, 0, 0, 0], [-1, 1, 0, 1], [0, -1, 1, 2]]
      for dx, dy, dz, idx in transition:
        s = px + py + pz
        dp[px][py][pz] += (get_expectation([px-dx, py-dy, pz-dz]) + N/s)*[px, py, pz][idx]/s
      
      return dp[px][py][pz]


print(get_expectation(pos))

Fixed submission:
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])

Buggy python submission (verdict: Time Limit Exceeded):
N = int(input())
A = list(map(int, input().split()))
px = py = pz = 0
for a in A:
  if a == 1:
    px += 1
  elif a == 2:
    py += 1
  elif a == 3:
    pz += 1

K = sum(A)

is_not_added = [[[True]*(N+1) for y in range(N+1)] for x in range(N+1)]
xyz = [[[0, 0, 0]]]
is_not_added[0][0][0] = False

transition = [[1, 0, 0, 0], [-1, 1, 0, 1], [0, -1, 1, 2]]
for k in range(1, K+1):
  _xyz = []
  for x, y, z in xyz[k-1]:
    for dx, dy, dz, _ in transition:
      if 0<=x+dx<=N and 0<=y+dy<=N and 0<=z+dz<=N \
          and is_not_added[x+dx][y+dy][z+dz]:
        is_not_added[x+dx][y+dy][z+dz] = False
        _xyz.append([x+dx, y+dy, z+dz])
  
  xyz.append(_xyz)


dp = [[[0]*(N+1) for y in range(N+1)] for x in range(N+1)]
for k in range(1, K+1):
  for x, y, z in xyz[k]:
    for dx, dy, dz, idx in transition:
      if 0<=x-dx<=N and 0<=y-dy<=N and 0<=z-dz<=N:
        w = [x, y, z][idx]
        s = x + y + z
        dp[x][y][z] += (dp[x-dx][y-dy][z-dz] + N/s) * w/s

print(dp[px][py][pz])

Fixed submission:
N = int(input())
A = input().split()
px = py = pz = 0
for a in A:
  if a == '1':
    px += 1
  elif a == '2':
    py += 1
  elif a == '3':
    pz += 1

nx = px + py + pz
ny = py + pz
nz = pz

dp = [[[0]*(nx+1) for y in range(ny+1)] for z in range(nz+1)]

for z in range(nz+1):
  for y in range(ny+1-z):
    for x in range(nx+1-z-y):
      if x == y == z == 0:
        continue
      
      else:
        s = x + y + z
        t = N/s
        if x > 0:
          t += dp[z][y][x-1]*x/s
        if y > 0:
          t += dp[z][y-1][x+1]*y/s
        if z > 0:
          t += dp[z-1][y+1][x]*z/s
        
        dp[z][y][x] += t


print(dp[pz][py][px])