Problem:
 Score : 100 points    Problem Statement There are N dishes, numbered 1, 2, \ldots, N .
Initially, for each i ( 1 \leq i \leq N ), Dish i has a_i ( 1 \leq a_i \leq 3 ) pieces of sushi on it.  Taro will perform the following operation repeatedly until all the pieces of sushi are eaten:   Roll a die that shows the numbers 1, 2, \ldots, N with equal probabilities, and let i be the outcome. If there are some pieces of sushi on Dish i , eat one of them; if there is none, do nothing.   Find the expected number of times the operation is performed before all the pieces of sushi are eaten.      Constraints  All values in input are integers.  1 \leq N \leq 300  1 \leq a_i \leq 3         Input Input is given from Standard Input in the following format:  N  a_1  a_2  \ldots  a_N       Output Print the expected number of times the operation is performed before all the pieces of sushi are eaten.
The output is considered correct when the relative difference is not greater than 10^{-9} .        Sample Input 1 3
1 1 1      Sample Output 1 5.5  The expected number of operations before the first piece of sushi is eaten, is 1 .
After that, the expected number of operations before the second sushi is eaten, is 1.5 .
After that, the expected number of operations before the third sushi is eaten, is 3 .
Thus, the expected total number of operations is 1 + 1.5 + 3 = 5.5 .       Sample Input 2 1
3      Sample Output 2 3  Outputs such as 3.00 , 3.000000003 and 2.999999997 will also be accepted.       Sample Input 3 2
1 2      Sample Output 3 4.5       Sample Input 4 10
1 3 2 3 3 2 3 2 1 3      Sample Output 4 54.48064457488221

Solution (Python):
import numpy as np
from numba import njit


@njit('(f8[:,:,:],i4,i4,i4,i4,)')
def solve(dp, N, n1, n2, n3):

    if dp[n1][n2][n3] >= 0.0:
        return dp[n1][n2][n3]
    if n1 == n2 == n3 == 0:
        return 0.0

    cnt = n1 + n2 + n3
    ret = 0
    if n3 >= 1:
        ret += solve(dp, N, n1, n2 + 1, n3 - 1) * n3
    if n2 >= 1:
        ret += solve(dp, N, n1 + 1, n2 - 1, n3) * n2
    if n1 >= 1:
        ret += solve(dp, N, n1 - 1, n2, n3) * n1

    ret = (ret + N) / cnt
    dp[n1][n2][n3] = ret
    return ret


def main():
    N = int(input())
    A = list(map(int, input().split()))

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

    a1, a2, a3 = 0, 0, 0
    for i in range(N):
        if A[i] == 1:
            a1 += 1
        elif A[i] == 2:
            a2 += 1
        elif A[i] == 3:
            a3 += 1

    print(solve(dp, N, a1, a2, a3))


if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np
from numba import njit


@njit(cache=True)
def solve(dp, N, n1, n2, n3):

    if dp[n1][n2][n3] >= 0.0:
        return dp[n1][n2][n3]
    if n1 == n2 == n3 == 0:
        return 0.0

    cnt = n1 + n2 + n3
    ret = 0
    if n3 >= 1:
        ret += solve(dp, N, n1, n2 + 1, n3 - 1) * n3
    if n2 >= 1:
        ret += solve(dp, N, n1 + 1, n2 - 1, n3) * n2
    if n1 >= 1:
        ret += solve(dp, N, n1 - 1, n2, n3) * n1

    ret = (ret + N) / cnt
    dp[n1][n2][n3] = ret
    return ret


def main():
    N = int(input())
    A = list(map(int, input().split()))

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

    a1, a2, a3 = 0, 0, 0
    for i in range(N):
        if A[i] == 1:
            a1 += 1
        elif A[i] == 2:
            a2 += 1
        elif A[i] == 3:
            a3 += 1

    print(solve(dp, N, a1, a2, a3))


if __name__ == '__main__':
    main()

Solution (Python):
n=int(input())
a=list(map(int,input().split()))
num1,num2,num3=0,0,0
for i in range(n):
    if a[i]==1:
        num1+=1
    elif a[i]==2:
        num2+=1
    else:
        num3+=1
dp=[[[-1.0 for k in range(n+1)] for i in range(n+1)] for j in range(n+1)]
def f(n1,n2,n3):
    if dp[n1][n2][n3]>=0.0:
        return dp[n1][n2][n3]
    if n1==n2==n3==0:
        return 0.0
    ret=n
    if n1>0:
        ret+=f(n1-1,n2,n3)*n1
    if n2>0:
        ret+=f(n1+1,n2-1,n3)*n2
    if n3>0:
        ret+=f(n1,n2+1,n3-1)*n3
    ret*=1.0/(n1+n2+n3)
    dp[n1][n2][n3]=ret
    return dp[n1][n2][n3]
print(f(num1,num2,num3))

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
from bisect import bisect_right, bisect_left
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, gamma, log
from operator import mul
from functools import reduce
from copy import deepcopy

sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.readline().split()))
def FI(): return list(map(float, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): pass
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 10 ** 9 + 7


n = I()
A = LI()
memo  = [[[-1.0] * (n+1) for _ in range(n+1)] for _ in range(n+1)]
memo[0][0][0] = 0
def f(i, j, k):
    if memo[i][j][k] != -1:
        return memo[i][j][k]
    s = i + j + k
    ret = n - s
    if i:
        ret += i * (f(i - 1, j, k) + 1)
    if j:
        ret += j * (f(i + 1, j - 1, k) + 1)
    if k:
        ret += k * (f(i, j + 1, k - 1) + 1)
    memo[i][j][k] = ret / s
    return memo[i][j][k]

one, two, three = 0, 0, 0
for i in range(n):
    if A[i] == 1:
        one += 1
    elif A[i] == 2:
        two += 1
    else:
        three += 1

print(f(one, two, three))

Solution (Python):
from collections import defaultdict

N = int(input())
A = list(map(int,input().split()))
dic = defaultdict(int)
for a in A:
  dic[a] += 1
a = dic[1]
b = dic[2]
c = dic[3]

dp = [[[0 for k in range(N+1)] for j in range(N+1)] for i in range(N+1)]
dp[a+b][a][0] = 1
for i in range(a+b,N+1):
  for j in range(a,i+1):
    for k in range(j+1):
      p,q,r = j-k,i-j,N-i
      if p:
        dp[i][j][k+1] += dp[i][j][k] * p/(p+q+r)
      if q:
        dp[i][j+1][k] += dp[i][j][k] * q/(p+q+r)
      if r:
        dp[i+1][j][k] += dp[i][j][k] * r/(p+q+r)

ans = 0
for i in range(a+b,N+1):
  for j in range(a,i+1):
    for k in range(j+1):
      p,q,r = j-k,i-j,N-i
      if p:
        ans += dp[i][j][k]*N/(p+q+r)*p/(p+q+r)
      if q:
        ans += dp[i][j][k]*N/(p+q+r)*q/(p+q+r)
      if r:
        ans += dp[i][j][k]*N/(p+q+r)*r/(p+q+r)
        
print(ans)

Solution (Python):






import sys
input = sys.stdin.readline
sys.setrecursionlimit(500000)

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]


def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()


def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count





N = I()
A = IL()
dic = DD(int)
for a in A:
  dic[a] += 1
a = dic[1]
b = dic[2]
c = dic[3]

dp = [[[0 for k in range(N+1)] for j in range(N+1)] for i in range(N+1)]
dp[a+b][a][0] = 1
for i in range(a+b,N+1):
  for j in range(a,i+1):
    for k in range(j+1):
      p,q,r = j-k,i-j,N-i
      if p:
        dp[i][j][k+1] += dp[i][j][k] * p/(p+q+r)
      if q:
        dp[i][j+1][k] += dp[i][j][k] * q/(p+q+r)
      if r:
        dp[i+1][j][k] += dp[i][j][k] * r/(p+q+r)

ans = 0
for i in range(a+b,N+1):
  for j in range(a,i+1):
    for k in range(j+1):
      p,q,r = j-k,i-j,N-i
      if p:
        ans += dp[i][j][k]*N/(p+q+r)*p/(p+q+r)
      if q:
        ans += dp[i][j][k]*N/(p+q+r)*q/(p+q+r)
      if r:
        ans += dp[i][j][k]*N/(p+q+r)*r/(p+q+r)
        
print(ans)

Solution (Python):
n=int(input())
a=list(map(int,input().split()))
num1,num2,num3=0,0,0
for i in range(n):
    if a[i]==1:
        num1+=1
    elif a[i]==2:
        num2+=1
    else:
        num3+=1
dp=[[[-1.0 for k in range(n+1)] for i in range(n+1)] for j in range(n+1)]
def f(n1,n2,n3):
    if dp[n1][n2][n3]>=0.0:
        return dp[n1][n2][n3]
    if n1==n2==n3==0:
        return 0.0
    ret=n
    if n1>0:
        ret+=f(n1-1,n2,n3)*n1
    if n2>0:
        ret+=f(n1+1,n2-1,n3)*n2
    if n3>0:
        ret+=f(n1,n2+1,n3-1)*n3
    ret*=1.0/(n1+n2+n3)
    dp[n1][n2][n3]=ret
    return dp[n1][n2][n3]
print(f(num1,num2,num3))

Solution (Python):
import sys
readline = sys.stdin.readline




N = int(readline())
A = list(map(int,readline().split()))

cnt = [0] * 4

for a in A:
  cnt[a] += 1

dp = [[[0 for k in range(N + 2)] for j in range(N + 2)] for i in range(N + 2)]

for i in range(N + 1):
  for j in range(N + 1):
    for k in range(N + 1):
      nonzero = i + j + k
      if nonzero == 0:
        continue
      dp[i][j][k] = N / nonzero
      if i - 1 >= 0:
        dp[i][j][k] += dp[i - 1][j + 1][k] * i / nonzero
      if j - 1 >= 0:
        dp[i][j][k] += dp[i][j - 1][k + 1] * j / nonzero
      if k - 1 >= 0:
        dp[i][j][k] += dp[i][j][k - 1] * k / nonzero

print(dp[cnt[3]][cnt[2]][cnt[1]])

Solution (Python):
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    
    
    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()

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)

N = int(input())
a = list(map(int, input().split()))

one = a.count(1)
two = a.count(2)
three = a.count(3)

dp = [[[0.0] * (N + 1) for _ in range((N + 1))] for _ in range((N + 1))]

max_one = one + two + three
max_two = two + three
max_three = three

for k in range(max_three + 1):
    for j in range(max_two + 1 - k):
        for i in range(max_one + 1 - j - k):
            if i + j + k == 0:
                continue
            else:
                r = N

            if i > 0:
                r += dp[i - 1][j][k] * i
            if j > 0:
                r += dp[i + 1][j - 1][k] * j
            if k > 0:
                r += dp[i][j + 1][k - 1] * k

            dp[i][j][k] = r / (i + j + k)

print(dp[one][two][three])

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)

N = int(input())
a = list(map(int, input().split()))

one = a.count(1)
two = a.count(2)
three = a.count(3)

dp = [[[-1] * 301 for _ in range(301)] for _ in range(301)]

def rec(i, j, k):
    if dp[i][j][k] >= 0:
        return dp[i][j][k]

    if i == 0 and j == 0 and k == 0:
        return 0.0

    res = 0.0

    if i > 0:
        res += rec(i - 1, j, k) * i
    if j > 0:
        res += rec(i + 1, j - 1, k) * j
    if k > 0:
        res += rec(i, j + 1, k - 1) * k
        
    res += N

    res /= i + j + k

    dp[i][j][k] = res
    return res

res = rec(one, two, three)

print(res)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)

N = int(input())
a = list(map(int, input().split()))

one = a.count(1)
two = a.count(2)
three = a.count(3)

dp = [[[-1] * 301 for _ in range(301)] for _ in range(301)]

def dfs(i, j, k):
    if dp[i][j][k] >= 0:
        return dp[i][j][k]

    if i == 0 and j == 0 and k == 0:
        return 0.0

    res = 0.0

    if i > 0:
        res += dfs(i - 1, j, k) * i
    if j > 0:
        res += dfs(i + 1, j - 1, k) * j
    if k > 0:
        res += dfs(i, j + 1, k - 1) * k
        
    res += N

    res /= i + j + k

    dp[i][j][k] = res
    return res

res = dfs(one, two, three)

print(res)

Solution (Python):
import sys
sys.setrecursionlimit(10**8)
def calc(N,i,j,k):
  if dp[i][j][k]>=0:
    return dp[i][j][k]
  ans = 1
  if k-1>=0 and j+1<=N:
    ans += calc(N,i,j+1,k-1)*k/N
  if j-1>=0 and i+1<=N:
    ans += calc(N,i+1,j-1,k)*j/N
  if i-1>=0:
    ans += calc(N,i-1,j,k)*i/N
  ans *= N/(i+j+k)
  dp[i][j][k] = ans
  return ans
from collections import Counter
N = int(input())
A = list(map(int, input().split()))
c = Counter(A)
dp = [[[-1]*(N+1) for _ in range(N+1)] for _ in range(N+1)]
dp[0][0][0] = 0
ans = calc(N,c[1],c[2],c[3])
print(ans)

Solution (Python):
n=int(input())
arr=list(map(int,input().split()))
one=two=three=0
for i in arr:
    if i==1:
        one+=1
    elif i==2:
        two+=1
    elif i==3:
        three+=1
max_one=one+two+three
max_two=two+three
max_three=three
dp=[[[0.0]*(n+1) for i in range(n+1)] for j in range(n+1)]
for k in range(max_three+1):
    for j in range(max_two+1-k):
        for i in range(max_one+1-j-k):
            if i+j+k==0:
                continue
            else:
                dp[i][j][k]=n/(i+j+k)
            if i>0:
                dp[i][j][k] += (dp[i-1][j][k]*i)/(i+j+k)
            if j>0:
                dp[i][j][k] += (dp[i+1][j-1][k]*j)/(i+j+k)
            if k>0:
                dp[i][j][k] += (dp[i][j+1][k-1]*k)/(i+j+k)
print(dp[i][j][k])

Solution (Python):
import collections

n=int(input())
arr=list(map(int,input().split()))
dp=[[[0]*(n+1) for _ in range(n+1)] for _ in range(n+1)]
cnt=collections.Counter(arr)
for total in range(1,n+1):
  for three in range(total+1):
    for two in range(total-three+1):
      one=total-two-three
      dp[one][two][three]+=n/total
      if one!=0:
        dp[one][two][three]+=dp[one-1][two][three]*one/total
      if two!=0:
        dp[one][two][three]+=dp[one+1][two-1][three]*two/total
      if three!=0:
        dp[one][two][three]+=dp[one][two+1][three-1]*three/total
print(dp[cnt[1]][cnt[2]][cnt[3]])

Solution (Python):
n = int(input())
A = list(map(int,input().split()))
one = A.count(1)
two = A.count(2)
three = A.count(3)

dp = [[[0]*(n+1) for _ in range(n+1)] for _ in range(n+1)]
for sushi in range(1,n+1):
    for k in range(sushi+1):
        for j in range(sushi-k+1):
            i = sushi-k-j
            if i != 0:
                dp[i][j][k] += dp[i - 1][j][k] * i / sushi
            if j != 0:
                dp[i][j][k] += dp[i + 1][j - 1][k] * j / sushi
            if k != 0:
                dp[i][j][k] += dp[i][j + 1][k - 1] * k / sushi
            dp[i][j][k] += n / sushi
print(dp[one][two][three])

Solution (Python):
import collections

n = int(input())
lst = list(map(int, input().split()))

c = collections.Counter(lst)
p, q, r = c[1], c[2], c[3]

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

for k in range(0, r + 1):
    for j in range(0, q + r  + 1 - k):
        for i in range(0, p + q + r + 1 - k - j):
            if i == 0 and j == 0 and k == 0:
                continue
            E = n
            if i > 0:
                E += dp[i-1][j][k] * i
            if j > 0:
                E += dp[i+1][j-1][k] * j
            if k  > 0:
                E += dp[i][j+1][k-1] * k
            dp[i][j][k] = E / (i + j + k)
print(dp[p][q][r])

Solution (Python):
from collections import Counter

def main():
    N = int(input())
    a = Counter(map(int, input().split()))
    a1, a2, a3 = a.get(1, 0), a.get(2, 0), a.get(3, 0)
    a123 = a1 + a2 + a3 + 1
    u = [[0.0] * (a123 - c2) for c2 in range(a123)]
    v = [[0.0] * (a123 - c2) for c2 in range(a123)]
    t1 = 0.0
    v2 = v[0]
    for c1 in range(1, a123):
        v2[c1] = t1 = N / c1 + t1
    for c2 in range(1, a123):
        v2 = v[c2]
        for c1, t2 in zip(range(301 - c2), v[c2 - 1][1:]):
            v2[c1] = t1 = (N + c1 * t1 + c2 * t2) / (c1 + c2)

    for c3 in range(1, a3 + 1):
        u, v = v, u
        v2 = v[0]
        for c1, t3 in zip(range(a123 - c3), u[1]):
            v2[c1] = t1 = (N + c1 * t1 + c3 * t3) / (c1 + c3)
        for c2 in range(1, a123 - c3):
            v2 = v[c2]
            c23 = c2 + c3
            for c1, t2, t3 in zip(range(a123 - c3 - c2), v[c2 - 1][1:], u[c2 + 1]):
                v2[c1] = t1 = (N + c1 * t1 + c2 * t2 + c3 * t3) / (c1 + c23)

    print(v[a2][a1])
main()

Solution (Python):

def resolve():
    N = int(input())
    A = list(map(int, input().split()))
    one = A.count(1)
    two = A.count(2)
    three = A.count(3)
    
    
    dp = [[[0.0] * (N+1) for _ in range(N+1)] for _ in range(N+1)]
    for c3 in range(N+1):
        for c2 in range(N+1):
            for c1 in range(N+1):
                if (c3 + c2 + c1) <= N:
                    
                    nozero = float(c3+c2+c1)
                    if nozero == 0.0:
                        continue
                    dp[c3][c2][c1] = float(N) / nozero
                    
                    
                    if c3-1 >= 0:
                        dp[c3][c2][c1] += dp[c3 - 1][c2+1][c1] * float(c3) / nozero
                    if c2-1 >= 0:
                        dp[c3][c2][c1] += dp[c3][c2 - 1][c1+1] * float(c2) / nozero
                    if c1-1 >= 0:
                        dp[c3][c2][c1] += dp[c3][c2][c1 - 1] * float(c1) / nozero

    print(dp[three][two][one])


if __name__ == "__main__":
    resolve()

Solution (Python):

def resolve():
    N = int(input())
    A = list(map(int, input().split()))
    one = A.count(1)
    two = A.count(2)
    three = A.count(3)

    dp = [[[0.0] * (N+1) for _ in range(N+1)] for _ in range(N+1)]
    for sushi in range(1, N + 1):
        for k in range(sushi + 1):
            for j in range(0, sushi - k + 1):
                i = sushi - k - j
                if i != 0:
                    dp[i][j][k] += dp[i - 1][j][k] * i / sushi
                if j != 0:
                    dp[i][j][k] += dp[i + 1][j - 1][k] * j / sushi
                if k != 0:
                    dp[i][j][k] += dp[i][j + 1][k - 1] * k / sushi
                dp[i][j][k] += N / sushi
    print(dp[one][two][three])


if __name__ == "__main__":
    resolve()

Solution (Python):
N = int(input())
a_list = list(map(int, input().split()))

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

cnt = {1:0, 2:0, 3:0}
for a in a_list:
    cnt[a] += 1

for i in range(N+1):
    for j in range(N+1):
        for k in range(N+1):
            if i==0 and j==0 and k==0:
                continue
            dp[i][j][k] += (N/(i+j+k))
            if k >= 1:
                dp[i][j][k] += dp[i][j][k-1] * (k/(i+j+k))
            if j >= 1:
                dp[i][j][k] += dp[i][j-1][k+1] * (j/(i+j+k))
            if i >= 1:
                dp[i][j][k] += dp[i-1][j+1][k] * (i/(i+j+k))
print(dp[cnt[3]][cnt[2]][cnt[1]])

Solution (Python):
import sys
def rec(N,dp,i,j,k):
    if dp[i][j][k]>=0.0:return dp[i][j][k]
    if all([i==0,j==0,k==0]):return 0.0
    res=0.0
    if i>0:res+=rec(N,dp,i-1,j,k)*i
    if j>0:res+=rec(N,dp,i+1,j-1,k)*j
    if k>0:res+=rec(N,dp,i,j+1,k-1)*k
    res=(res+N)*1.0/(i+j+k)
    dp[i][j][k]=res
    return dp[i][j][k]
def main():
    N=int(sys.stdin.readline())
    a=tuple(map(int,sys.stdin.readline().split()))
    n1,n2,n3=a.count(1),a.count(2),a.count(3)
    dp=[[[-1.0 for _ in range(n3+1)] for _ in range(n2+n3+1)] for _ in range(n1+n2+n3+1)]
    print(rec(N,dp,n1,n2,n3))
if __name__=='__main__':main()

Solution (Python):
import sys
sys.setrecursionlimit(10**8)
stdin = sys.stdin


def ni(): return int(ns())


def na(): return list(map(int, stdin.readline().split()))


def naa(N): return [na() for _ in range(N)]


def ns(): return stdin.readline().rstrip()  


N = ni()
a_array = na()
cnt = [0] * 3
for a in a_array:
    cnt[a-1] += 1

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

for sushi in range(1, N+1):
    for k in range(sushi+1):
        for j in range(0, sushi-k+1):
            i = sushi-k-j
            if i != 0:
                dp[i][j][k] += dp[i-1][j][k] * i / sushi
            if j != 0:
                dp[i][j][k] += dp[i+1][j-1][k] * j / sushi
            if k != 0:
                dp[i][j][k] += dp[i][j+1][k-1] * k / sushi
            dp[i][j][k] += N / sushi

print(dp[cnt[0]][cnt[1]][cnt[2]])

Solution (Python):
import sys
from collections import Counter
sys.setrecursionlimit(10000)

N = int(input())
a = list(map(int, input().split()))
dp = [[[-1 for i in range(N+1)] for j in range(N+1)] for k in range(N+1)]

def rec(i, j, k):
    if dp[i][j][k] >= 0:
        return dp[i][j][k]
    if i == 0 and j == 0 and k == 0:
        return 0.0
    res = 0.0
    if i > 0:
        res += rec(i-1, j, k) * i
    if j > 0:
        res += rec(i+1, j-1, k) * j
    if k > 0:
        res += rec(i, j+1, k-1) * k
    res += N
    res *= 1.0/(i+j+k)

    dp[i][j][k] = res
    return res

c = Counter(a)
one, two, three = c[1], c[2], c[3]
print(rec(one, two, three))

Solution (Python):


import sys

sys.setrecursionlimit(10**6)
input = sys.stdin.buffer.readline
INF = 10 ** 9 + 1  


def debug(*x):
    print(*x)


def solve(N, AS):
    from collections import Counter
    count = Counter(AS)
    expect = [[[-1] * (N + 1) for i in range(N + 1)] for j in range(N + 1)]
    expect[0][0][0] = 0

    def f(a, b, c):
        p = N
        if a > 0:
            v = expect[a - 1][b][c]
            if v == -1:
                v = f(a - 1, b, c)
            p += v * a
        if b > 0:
            v = expect[a + 1][b - 1][c]
            if v == -1:
                v = f(a + 1, b - 1, c)
            p += v * b
        if c > 0:
            v = expect[a][b + 1][c - 1]
            if v == -1:
                v = f(a, b + 1, c - 1)
            p += v * c
        p /= (a + b + c)

        
        expect[a][b][c] = p
        return p

    return f(count[1], count[2], count[3])


def main():
    N = int(input())
    AS = list(map(int, input().split()))
    print(solve(N, AS))


T1 = """
3
1 1 1
"""


def test_T1():
    """
    >>> as_input(T1)
    >>> main()
    5.5
    """


T2 = """
1
3
"""


def test_T2():
    """
    >>> as_input(T2)
    >>> main()
    3.0
    """


T3 = """
2
1 2
"""


def test_T3():
    """
    >>> as_input(T3)
    >>> main()
    4.5
    """

















def _test():
    import doctest
    doctest.testmod()


def as_input(s):
    "use in test, use given string as input file"
    import io
    global read, input
    f = io.StringIO(s.strip())
    input = f.readline
    read = f.read


USE_NUMBA = False
if (USE_NUMBA and sys.argv[-1] == 'ONLINE_JUDGE') or sys.argv[-1] == '-c':
    print("compiling")
    from numba.pycc import CC
    cc = CC('my_module')
    cc.export('solve', solve.__doc__.strip().split()[0])(solve)
    cc.compile()
    exit()
else:
    input = sys.stdin.buffer.readline
    read = sys.stdin.buffer.read

    if (USE_NUMBA and sys.argv[-1] != '-p') or sys.argv[-1] == "--numba":
        
        
        from my_module import solve  
    elif sys.argv[-1] == "-t":
        _test()
        sys.exit()
    elif sys.argv[-1] != '-p' and len(sys.argv) == 2:
        
        input_as_file = open(sys.argv[1])
        input = input_as_file.buffer.readline
        read = input_as_file.buffer.read

    main()

Solution (Python):
import sys
from collections import Counter
sys.setrecursionlimit(10000)

N = int(input())
a = list(map(int, input().split()))
dp = [[[-1 for i in range(N+1)] for j in range(N+1)] for k in range(N+1)]

def rec(i, j, k):
    if dp[i][j][k] >= 0:
        return dp[i][j][k]
    if i == 0 and j == 0 and k == 0:
        return 0.0
    res = 0.0
    if i > 0:
        res += rec(i-1, j, k) * i
    if j > 0:
        res += rec(i+1, j-1, k) * j
    if k > 0:
        res += rec(i, j+1, k-1) * k
    res += N
    res *= 1.0/(i+j+k)

    dp[i][j][k] = res
    return res

c = Counter(a)
one, two, three = c[1], c[2], c[3]
print(rec(one, two, three))

Solution (Python):
import collections

n = int(input())
lst = list(map(int, input().split()))

c = collections.Counter(lst)
p, q, r = c[1], c[2], c[3]

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

for k in range(0, r + 1):
    for j in range(0, q + r  + 1 - k):
        for i in range(0, p + q + r + 1 - k - j):
            if i == 0 and j == 0 and k == 0:
                continue
            E = n
            if i > 0:
                E += dp[i-1][j][k] * i
            if j > 0:
                E += dp[i+1][j-1][k] * j
            if k  > 0:
                E += dp[i][j+1][k-1] * k
            dp[i][j][k] = E / (i + j + k)
print(dp[p][q][r])

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n = int(input())
a = list(map(int, input().split()))

from collections import Counter, defaultdict
c = Counter(a)
n1, n2, n3 = c[1], c[2], c[3]
n0 = c[0]
def glist(*l, fill=None):
    if len(l)==1:
        return [fill]*l[0]
    ll = l[1:]
    return [glist(*ll) for _ in range(l[0])]

count = 0
S = sum(a)
sub = glist(n1+n2+n3+1,n2+n3+1,n3+1)
sub[0][0][0] = 0


for i3 in range(n3+1):
    for i2 in range(n2+n3+1-i3):
        for i1 in range(n1+n2+n3 - (i2+i3) + 1):
            if i1==i2==i3==0 or i1+i2+i3>n:
                continue

            ans = n
            p = 0
            if i1>0:
                ans += i1 * sub[i1-1][i2][i3]
                p += i1
            if i2>0:
                ans += i2 * sub[i1+1][i2-1][i3]
                p += i2
            if i3>0:
                ans += i3 * sub[i1][i2+1][i3-1]
                p += i3
            ans = ans/p
            sub[i1][i2][i3] = ans
print(sub[n1][n2][n3])

Solution (Python):
N = int(input())
a = list(map(int, input().split()))

one, two, three = 0, 0, 0
for ai in a:
    if ai == 1: one += 1
    if ai == 2: two += 1
    if ai == 3: three += 1

u = N + 1
dp = [[[-1]* u for _ in range(u)] for _ in range(u)]
dp[0][0][0] = 0

def rec(i, j, k):
    if dp[i][j][k] >= 0:
        return dp[i][j][k]
    res = N
    if i > 0: res += rec(i-1, j, k) * i
    if j > 0: res += rec(i+1, j-1, k) * j
    if k > 0: res += rec(i, j+1, k-1) * k
    res /= (i + j + k)
    dp[i][j][k] = res
    return res

print(rec(one, two, three))

Solution (Python):
N = int(input())

dp = [[[0.0]*(N+1) for _ in range(N+1)] for _ in range(N+1)]

dic = {}
dic[0] = 0
dic[1] = 0
dic[2] = 0
dic[3] = 0
for a in map(int, input().split(" ")):
    dic[a] += 1


for k in range(dic[3]+1):
    for j in range(dic[2]+1+dic[3]-k):
        for i in range(dic[1]+1+dic[2]-j+dic[3]-k):
            if i == 0 and j == 0 and k == 0:
                continue
            res = 0.0
            if i > 0:
                res += dp[i-1][j][k] * i
            if j > 0:
                res += dp[i+1][j-1][k] * j
            if k > 0:
                res += dp[i][j+1][k-1] * k
            res += N
            res *= (1.0 / (i+j+k))
            dp[i][j][k] = res

print(dp[dic[1]][dic[2]][dic[3]])

Solution (Python):
from collections import Counter

def solve():
    N = int(input())
    As = list(map(int, input().split()))

    cnt = Counter(As)

    dp = [[[0.0]*(N+1) for j in range(N+1)] for i in range(N+1)]
    for k in range(N+1):
        for j in range(N+1):
            if j+k > N:
                break
            for i in range(N+1):
                if i+j+k > N:
                    break
                expe = 0.0
                if i > 0:
                    expe += i/N * dp[i-1][j][k]
                if j > 0:
                    expe += j/N * dp[i+1][j-1][k]
                if k > 0:
                    expe += k/N * dp[i][j+1][k-1]
                if i+j+k:
                    expe = (expe+1)*N/(i+j+k)
                dp[i][j][k] = expe

    ans = dp[cnt[1]][cnt[2]][cnt[3]]
    print(ans)


solve()

Solution (Python):
import sys
input = sys.stdin.readline
n=int(input())
a=list(map(int,input().split()))

A=a.count(1)
B=a.count(2)
C=a.count(3)

dp=[[[0.0]*(n+1) for _ in range(n+1)] for _ in range(n+1)]
def c():
    for k in range(C+1):
        for j in range(B+C+1-k):
            for i in range(n+1-k-j):
                if i==0 and j==0 and k==0:
                    continue
                res=1.0*n
                if i>0:
                    res+=i*dp[i-1][j][k]
                if j>0:
                    res+=j*dp[i+1][j-1][k]
                if k>0:
                    res+=k*dp[i][j+1][k-1]
                dp[i][j][k]=res/(i+j+k)
    print(dp[A][B][C])
c()