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):
from sys import setrecursionlimit
setrecursionlimit(10**8)

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

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

    def search(i, j, k):
        num = n
        if i > 0:
            if dp[i-1][j][k] == -1:
                num += search(i-1, j, k) * i
            else:
                num += dp[i-1][j][k] * i
        if j > 0:
            if dp[i+1][j-1][k] == -1:
                num += search(i+1, j-1, k) * j
            else:
                num += dp[i+1][j-1][k] * j
        if k > 0:
            if dp[i][j+1][k-1] == -1:
                num += search(i, j+1, k-1) * k
            else:
                num += dp[i][j+1][k-1] * k
        num /= i + j + k

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

    print(search(a.count(1), a.count(2), a.count(3)))

main()

Solution (Python):
from collections import Counter
import sys


def main():
    sys.setrecursionlimit(10**6)
    readline = sys.stdin.buffer.readline
    N = int(readline())
    A = list(map(int, readline().split()))
    C = Counter(A)

    l, m, n = C[3], C[2], C[1]
    dp = [[[-1]*(N+1) for _ in range(N+1)] for _ in range(N+1)]
    dp[0][0][0] = 0

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


    print(exval(l, m, n))
    
if __name__ == "__main__":
    main()

Solution (Python):

import sys
sys.setrecursionlimit(10**6)
readline = sys.stdin.buffer.readline


def read_int():
    return int(readline())


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


def f(c1, c2, c3):
    if memo[c1][c2][c3] != -1:
        return memo[c1][c2][c3]
    c0 = N - (c1 + c2 + c3)
    p0 = c0 / N
    ip0 = 1 - p0
    ret = 1 / ip0
    if c1 > 0:
        p1 = c1 / N
        ret += f(c1-1, c2, c3) * p1 / ip0
    if c2 > 0:
        p2 = c2 / N
        ret += f(c1+1, c2-1, c3) * p2 / ip0
    if c3 > 0:
        p3 = c3 / N
        ret += f(c1, c2+1, c3-1) * p3 / ip0
    memo[c1][c2][c3] = ret
    return ret


N = read_int()
A = read_int_n()
memo = [[[-1] * (N + 1) for _ in range(N + 1)] for _ in range(N + 1)]
memo[0][0][0] = 0
print(f(A.count(1), A.count(2), A.count(3)))

Solution (Python):


import sys
sys.setrecursionlimit(10 ** 9)


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

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

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

def calc(a, b, c):
    if dp[a][b][c] >= 0:
        return dp[a][b][c]

    if a + b + c == 0:
        return 0

    t = N

    if a > 0:
        t += calc(a - 1, b + 1, c) * a
    if b > 0:
        t += calc(a, b - 1, c + 1) * b
    if c > 0:
        t += calc(a, b, c - 1) * c
        
    t /= a + b + c

    dp[a][b][c] = t
    return t

ans = calc(three, two, one)
print(ans)

Solution (Python):
from collections import Counter

N = int(input())
A = list(map(int, input().split()))
B = Counter(A)
DP = [[[0.0] * (B[3] + 1) for _ in range(N + 1)] for __ in range(N + 2)]

for k in range(B[3]+1):
    for j in range(B[2] + B[3] + 1 - k):
        for i in range(B[1] + B[2] + B[3] + 1 - k - j):
            if i == j == k == 0:
                continue
            tmp = float(N)
            if i > 0:
                tmp += i * DP[i - 1][j][k]
            if j > 0:
                tmp += j * DP[i + 1][j - 1][k]
            if k > 0:
                tmp += k * DP[i][j + 1][k - 1]
            DP[i][j][k] = tmp / (i + j + k)
print(DP[B[1]][B[2]][B[3]])

Solution (Python):
N = int(input())
a = list(map(float, input().split()))
dp=[[[0.0] *(N+2) for _ in range(N+2)] for _ in range(N+2)]
one,two,three=a.count(1),a.count(2),a.count(3)
for k in range(three+1):
    for j in range(two+three+1):
        for i in range(N+1):
            if k==0 and j==0 and i==0:
                continue
            num=N
            if i>0:
                num+=dp[i-1][j][k]*i
            if j>0:
                num+=dp[i+1][j-1][k]*j
            if k>0:
                num+=dp[i][j+1][k-1]*k
            dp[i][j][k]=num/(k+j+i)
print(dp[one][two][three])

Solution (Python):
N = int(input())
a = list(map(float, input().split()))
dp=[[[0.0] *(N+2) for _ in range(N+2)] for _ in range(N+2)]
one,two,three=a.count(1),a.count(2),a.count(3)
for k in range(three+1):
    for j in range(two+three+1-k):
        for i in range(N+1):
            if k==0 and j==0 and i==0:
                continue
            num=N
            if i>0:
                num+=dp[i-1][j][k]*i
            if j>0:
                num+=dp[i+1][j-1][k]*j
            if k>0:
                num+=dp[i][j+1][k-1]*k
            dp[i][j][k]=num/(k+j+i)
print(dp[one][two][three])

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

def main():
    n = int(input())
    alst = list(map(float, input().split()))
    dp = [[[0 for _ in range(n + 1)] for _ in range(n + 1)] for _ in range(n + 1)]
    one = alst.count(1)
    two = alst.count(2)
    three = alst.count(3)
    for k in range(three + 1):
        for j in range(two + three + 1):
            if k + j > n:
                continue
            for i in range(n + 1):
                if k + j + i > n:
                    continue
                ijk = i + j + k
                if ijk == 0:
                    continue
                tmp = n
                if i != 0:
                    tmp += dp[i - 1][j][k] * i
                if j != 0 and i != n:
                    tmp += dp[i + 1][j - 1][k] * j
                if k != 0 and j != n:
                    tmp += dp[i][j + 1][k - 1] * k
                dp[i][j][k] = tmp / ijk
    
    print(dp[one][two][three])


main()

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


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

for k in range(three+1):
  for j in range(two+three+1-k):
    for i in range(one+two+three+1-k-j):
      if i==0 and j==0 and k==0:
        continue
      tmp=n*1.0
      if i>0:
        tmp+=i*dp[i-1][j][k]
      if j>0:
        tmp+=j*dp[i+1][j-1][k]
      if k>0:
        tmp+=k*dp[i][j+1][k-1]
      dp[i][j][k]=tmp/(i+j+k)

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

Solution (Python):
n = int(input())
a = list(map(int,input().split()))
d = [0 for i in range(3)]
dp = [[[0 for i in range(n+1)] for j in range(n+1)] for k in range(n+1)] 
for i in a:
    d[i-1] +=1
for k in range(n+1):
    for j in range(n+1):
        for i in range(n+1):
            if i+j+k==0 or i+j+k>n:
                continue
            ans = n/(i+j+k)
            x = i/(i+j+k)
            if x:
                ans += x*dp[k][j][i-1]
            x = j/(i+j+k)
            if x:
                ans += x*dp[k][j-1][i+1]
            x = k/(i+j+k)
            if x:
                ans += x*dp[k-1][j+1][i]
            dp[k][j][i] = ans
print(dp[d[2]][d[1]][d[0]])

Solution (Python):


def main():
    import sys
    sys.setrecursionlimit(10**6)
    N = int(input())
    A = [int(i) for i in input().split()]
    dp = [0]*(N*((N+1)**2) + N*(N+1) + N + 1)
    c = [0]*4
    for a in A:
        c[a] += 1

    def rec(c1, c2, c3):
        if c1 == c2 == c3 == 0:
            return 0
        if dp[c1 * (N+1)**2 + c2 * (N+1) + c3] != 0:
            return dp[c1 * (N+1)**2 + c2 * (N+1) + c3]
        ret = N
        if c1 > 0:
            ret += rec(c1-1, c2, c3) * c1
        if c2 > 0:
            ret += rec(c1 + 1, c2-1, c3) * c2
        if c3 > 0:
            ret += rec(c1, c2 + 1, c3-1) * c3
        ret /= c1 + c2 + c3
        dp[c1 * (N+1)**2 + c2 * (N+1) + c3] = ret
        return ret

    print(rec(*c[1:]))


if __name__ == '__main__':
    main()

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)
n = int(input())
A = tuple(map(int, input().split()))
C = [0]*3
for a in A:
  C[a-1] += 1
dp = [[[-1]*(n+1) for i in range(n+1)] for j in range(n+1)]
dp[0][0][0] = 0
def f(i, j, k):
  if dp[i][j][k] != -1:
    return dp[i][j][k]
  res = 1
  if i:
    res += f(i-1, j, k)*(i/n)
  if j:
    res += f(i+1, j-1, k)*(j/n)
  if k:
    res += f(i, j+1, k-1)*(k/n)
  res /= (i+j+k)/n
  dp[i][j][k] = res
  return res
ans = f(C[0], C[1], C[2])
print(ans)

Solution (Python):


from __future__ import division, print_function
from fractions import Fraction
import sys
import os
from io import BytesIO, IOBase




from bisect import bisect_right as br
from bisect import bisect







BUFSIZE = 8192


class FastIO(IOBase):
    newlines = 0

    def __init__(self, file):
        self._fd = file.fileno()
        self.buffer = BytesIO()
        self.writable = "x" in file.mode or "r" not in file.mode
        self.write = self.buffer.write if self.writable else None

    def read(self):
        while True:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            if not b:
                break
            ptr = self.buffer.tell()
            self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr)
        self.newlines = 0
        return self.buffer.read()

    def readline(self):
        while self.newlines == 0:
            b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
            self.newlines = b.count(b"\n") + (not b)
            ptr = self.buffer.tell()
            self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr)
        self.newlines -= 1
        return self.buffer.readline()

    def flush(self):
        if self.writable:
            os.write(self._fd, self.buffer.getvalue())
            self.buffer.truncate(0), self.buffer.seek(0)


class IOWrapper(IOBase):
    def __init__(self, file):
        self.buffer = FastIO(file)
        self.flush = self.buffer.flush
        self.writable = self.buffer.writable
        self.write = lambda s: self.buffer.write(s.encode("ascii"))
        self.read = lambda: self.buffer.read().decode("ascii")
        self.readline = lambda: self.buffer.readline().decode("ascii")


def print(*args, **kwargs):
    """Prints the values to a stream, or to sys.stdout by default."""
    sep, file = kwargs.pop("sep", " "), kwargs.pop("file", sys.stdout)
    at_start = True
    for x in args:
        if not at_start:
            file.write(sep)
        file.write(str(x))
        at_start = False
    file.write(kwargs.pop("end", "\n"))
    if kwargs.pop("flush", False):
        file.flush()


if sys.version_info[0] < 3:
    sys.stdin, sys.stdout = FastIO(sys.stdin), FastIO(sys.stdout)
else:
    sys.stdin, sys.stdout = IOWrapper(sys.stdin), IOWrapper(sys.stdout)







mod = 1000000007


def inp(): return sys.stdin.readline().rstrip("\r\n")  


def out(var): sys.stdout.write(str(var))  


def lis(): return list(map(int, inp().split()))


def stringlis(): return list(map(str, inp().split()))


def sep(): return map(int, inp().split())


def strsep(): return map(str, inp().split())


def zerolist(n): return [0] * n


def nextline(): out("\n")  


def testcase(t):
    for p in range(t):
        solve()


from collections import defaultdict, deque


class Graph(object):
    def __init__(self, dic):
        self.dict = dic
        self.vertex = set()

    def addedge(self, a, b):
        if self.dict[a] == []:
            self.vertex.add(a)
        if b != None:
            self.dict[a].append(b)

    def con_com_util(self, i, unvisited, count, comp):
        for j in self.dict[i]:
            if j in unvisited:
                unvisited.remove(j)
                comp[count].append(j)
                comp = self.con_com_util(j, unvisited, count, comp)
        return comp

    def conn_components(self):
        unvisited = self.vertex.copy()
        count = 0
        comp = defaultdict(list)
        while unvisited:
            i = unvisited.pop()
            comp[count].append(i)
            comp = self.con_com_util(i, unvisited, count, comp)
            count += 1
        return (count, comp)


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

    return count



from sys import setrecursionlimit
setrecursionlimit(1000000)


def find(a0,a1,a2,a3):
    if a1+a2+a3==0:
        return 0
    if d[(a0,a1,a2,a3)]!=0:
        return d[(a0,a1,a2,a3)]
    if a1<0 or a2<0 or a3<0:
        return -1
    p0=a0/n
    p1=a1/n
    p2=a2/n
    p3=a3/n
    d[(a0, a1, a2, a3)]= p0 + p1*(find(a0+1,a1-1,a2,a3)+1) + p2*(find(a0,a1+1,a2-1,a3)+1) +p3*(find(a0,a1,a2+1,a3-1)+1)
    d[(a0, a1, a2, a3)]/=(1-p0)
    return d[(a0,a1,a2,a3)]





n=int(inp())
temp=lis()
ar=[0,0,0,0]
dp=[[[0 for _ in range(n+5)] for i in range(n+5)] for j in range(n+5)]
for i in temp:
    if i==1:
        ar[1]+=1
    elif i==2:
        ar[2]+=1
    else:
        ar[3]+=1
ma_3=ar[3]
ma_2=ar[2]+ar[3]
ma_1=ar[1]+ar[2]+ar[3]
for a3 in range(ma_3 +1):
    for a2 in range(ma_2 + 1 - a3):
        for a1 in range(ma_1+ 1 -a3-a2):
            if a1+a2+a3==0:
                continue
            a0=n-a1-a2-a3
            p0 = a0 / n
            p1 = a1 / n
            p2 = a2 / n
            p3 = a3 / n
            dp[a1][a2][a3] +=p0
            if a1!=0:
                dp[a1][a2][a3]+=(p1*(dp[a1-1][a2][a3]+1))
            if a2!=0:
                dp[a1][a2][a3]+=(p2*(dp[a1+1][a2-1][a3]+1))
            if a3!=0:
                dp[a1][a2][a3]+=(p3*(dp[a1][a2+1][a3-1]+1))
            dp[a1][a2][a3]/=(1-p0)
print(dp[ar[1]][ar[2]][ar[3]])

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)

ra = range
enu = enumerate


def exit(*argv, **kwarg):
    print(*argv, **kwarg)
    sys.exit()


def mina(*argv, sub=1): return list(map(lambda x: x - sub, argv))



def a_int(): return int(readline())


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W):
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


def read_col(H):
    '''H is number of rows
    A列、B列が与えられるようなとき
    ex1)A,B=read_col(H)    ex2) A,=read_col(H) 
    ret = []
    for _ in range(H):
        ret.append(list(map(int, readline().split())))
    return tuple(map(list, zip(*ret)))


def read_tuple(H):
    '''H is number of rows'''
    ret = []
    for _ in range(H):
        ret.append(tuple(map(int, readline().split())))
    return ret


MOD = 10**9 + 7
INF = 2**31  


from collections import defaultdict, Counter, deque
from operator import itemgetter, xor, add
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right  
from functools import reduce, lru_cache
from math import gcd


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


N = a_int()
A = ints()
cnt = Counter(A)
'''
dp[i,j,k]...すべての寿司がなくなるまでの操作の期待値(寿司が1つ乗っている皿の数がi,2つがj,3つがk)

dp[i,j,k]= 1(操作回数)+ 
            i/N*(1つの皿が一つ減った場合の操作の期待値) + 
            j/N*(2つの皿が1つ減った場合の操作の期待値) + 
            k/N*(3つの皿が1つ減った場合の操作の期待値)+
            (1-i-j-k)/N * (0つの皿を選んでしまったときの期待値)

dp[i,j,k]= 1+ i/N * dp[i-1,j,k]+ 
            j/N*dp[i+1,j-1,k] + 
            k/N*dp[i,j+1,k-1]+
            (N-i-j-k)/N*dp[i,j,k]

dp[i,j,k]= N/(i+j+k) *
            (1+ i/N * dp[i-1,j,k]+ 
            j/N*dp[i+1,j-1,k] + 
            k/N*dp[i,j+1,k-1])
'''

dp = np.full((N + 1, N + 1, N + 1), -1.0)
dp[0, 0, 0] = 0
from numba import njit



@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def dfs(i, j, k, N, dp):
    if i < 0 or j < 0 or k < 0 or i == j == k == 0:
        return 0
    if dp[i, j, k] != -1:
        return dp[i, j, k]
    ret = 1
    ret += i / N * dfs(i - 1, j, k, N, dp)
    ret += j / N * dfs(i + 1, j - 1, k, N, dp)
    ret += k / N * dfs(i, j + 1, k - 1, N, dp)
    if (i + j + k) == 0:
        dp[i, j, k] = ret
        return ret
    dp[i, j, k] = N / (i + j + k) * ret
    return dp[i, j, k]


print(dfs(cnt[1], cnt[2], cnt[3], N, dp))

Solution (Python):
n=int(input())
l=list(map(int,input().split()))
ones,twos,threes=0,0,0
for i in l:
    if i==1:
        ones+=1
    elif i==2:
        twos+=1
    else:
        threes+=1
max_ones=ones+twos+threes
max_twos=twos+threes
max_threes=threes
dp=[[[0.0]*(n+1)for i in range(n+1)]for j in range(n+1)]
for k in range(max_threes+1):
    for j in range(max_twos+1-k):
        for i in range(max_ones+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 sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)

ra = range
enu = enumerate


def exit(*argv, **kwarg):
    print(*argv, **kwarg)
    sys.exit()


def mina(*argv, sub=1): return list(map(lambda x: x - sub, argv))



def a_int(): return int(readline())


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W):
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


def read_col(H):
    '''H is number of rows
    A列、B列が与えられるようなとき
    ex1)A,B=read_col(H)    ex2) A,=read_col(H) 
    ret = []
    for _ in range(H):
        ret.append(list(map(int, readline().split())))
    return tuple(map(list, zip(*ret)))


def read_tuple(H):
    '''H is number of rows'''
    ret = []
    for _ in range(H):
        ret.append(tuple(map(int, readline().split())))
    return ret


MOD = 10**9 + 7
INF = 2**31  


from collections import defaultdict, Counter, deque
from operator import itemgetter, xor, add
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right  
from functools import reduce, lru_cache
from math import gcd


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


N = a_int()
A = ints()
cnt = Counter(A)
'''
dp[i,j,k]...すべての寿司がなくなるまでの操作の期待値(寿司が1つ乗っている皿の数がi,2つがj,3つがk)

dp[i,j,k]= 1(操作回数)+ 
            i/N*(1つの皿が一つ減った場合の操作の期待値) + 
            j/N*(2つの皿が1つ減った場合の操作の期待値) + 
            k/N*(3つの皿が1つ減った場合の操作の期待値)+
            (1-i-j-k)/N * (0つの皿を選んでしまったときの期待値)

dp[i,j,k]= 1+ i/N * dp[i-1,j,k]+ 
            j/N*dp[i+1,j-1,k] + 
            k/N*dp[i,j+1,k-1]+
            (N-i-j-k)/N*dp[i,j,k]

dp[i,j,k]= N/(i+j+k) *
            (1+ i/N * dp[i-1,j,k]+ 
            j/N*dp[i+1,j-1,k] + 
            k/N*dp[i,j+1,k-1])
'''

dp = np.full((N + 1, N + 1, N + 1), -1.0)
dp[0, 0, 0] = 0
from numba import njit



@njit('(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def dfs(i, j, k, N, dp):
    if i < 0 or j < 0 or k < 0 or i == j == k == 0:
        return 0
    if dp[i, j, k] != -1:
        return dp[i, j, k]
    ret = 1
    ret += i / N * dfs(i - 1, j, k, N, dp)
    ret += j / N * dfs(i + 1, j - 1, k, N, dp)
    ret += k / N * dfs(i, j + 1, k - 1, N, dp)
    if (i + j + k) == 0:
        dp[i, j, k] = ret
        return ret
    dp[i, j, k] = N / (i + j + k) * ret
    return dp[i, j, k]


print(dfs(cnt[1], cnt[2], cnt[3], N, dp))

Solution (Python):
import sys
from collections import Counter as cc
from numba import njit
import numpy as np

def input(): return sys.stdin.readline().rstrip()
def ii(): return int(input())
def mi(): return map(int, input().split())
def li(): return list(mi())


n = ii()
a = li()
c = cc(a)
N = sum(c.values())
cache = np.full((N+1, N+1, N+1), -1, dtype=float)
cache[0][0][0] = 0
@njit(cache=True)
def dfs(x, y, z, N, cache):
    if x+y+z==0:
        return 0
    if cache[x][y][z]>=0:
        return cache[x][y][z]
    tmp = 1
    if x > 0:
        tmp += dfs(x-1, y, z,N, cache)*x/N
    if y > 0:
        tmp += dfs(x+1, y-1, z, N, cache)*y/N
    if z > 0:
        tmp += dfs(x, y+1, z-1, N, cache)*z/N
    tmp *= N/(x+y+z)
    cache[x][y][z] = tmp
    return tmp

dfs(c[1], c[2], c[3], N, cache)
print(cache[c[1]][c[2]][c[3]])

Solution (Python):
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])

Solution (Python):
import sys
from collections import Counter
input = sys.stdin.readline


N = int(input())
memo = [[[-1]*(N+1) for _ in range(N+1)] for _ in range(N+1)]

def dp(a, b, c):
    if a == b == c == 0: return 0
    if memo[a][b][c] != -1:
        return memo[a][b][c]
    
    res = N
    if a > 0: res += a*dp(a-1, b, c)
    if b > 0: res += b*dp(a+1, b-1, c)
    if c > 0: res += c*dp(a, b+1, c-1)
    res /= a+b+c
    
    memo[a][b][c] = res
    return res

def main():
    A = list(map(int, input().split()))
    
    count = Counter(A)
    for i in range(1, 4):
        if not i in count: count[i] = 0
    x, y, z = count[1], count[2], count[3]

    
    ans = dp(x, y, z)
    
    print(ans)
        

    
if __name__ == "__main__":
    main()

Solution (Python):
N = int(input())
A = list(map(int, input().split()))
cnt1 = A.count(1)
cnt2 = A.count(2)
cnt3 = A.count(3)






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


def func(c1, c2, c3):
    if dp[c1][c2][c3] != -1:
        return dp[c1][c2][c3]

    prob0 = (N - c1 - c2 - c3) / N
    temp = 1 / (1 - prob0)
    if c1 > 0:
        temp += func(c1 - 1, c2, c3) * (c1 / N) / (1 - prob0)
    if c2 > 0:
        temp += func(c1 + 1, c2 - 1, c3) * (c2 / N) / (1 - prob0)
    if c3 > 0:
        temp += func(c1, c2 + 1, c3 - 1) * (c3 / N) / (1 - prob0)

    dp[c1][c2][c3] = temp
    return temp


print(func(cnt1, cnt2, cnt3))

Solution (Python):
N = int(input())
A = list(map(int, input().split()))
cnt1 = A.count(1)
cnt2 = A.count(2)
cnt3 = A.count(3)

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

def func(c1, c2, c3):
    if dp[c1][c2][c3] != -1:
        return dp[c1][c2][c3]

    prob0 = (N - c1 - c2 - c3) / N
    temp = 1 / (1 - prob0)
    if c1 > 0:
        temp += func(c1 - 1, c2, c3) * (c1 / N) / (1 - prob0)
    if c2 > 0:
        temp += func(c1 + 1, c2 - 1, c3) * (c2 / N) / (1 - prob0)
    if c3 > 0:
        temp += func(c1, c2 + 1, c3 - 1) * (c3 / N) / (1 - prob0)

    dp[c1][c2][c3] = temp
    return temp

print(func(cnt1, cnt2, cnt3))

Solution (Python):



N = int(input())
A = list(map(int, input().split()))
cnt1 = A.count(1)
cnt2 = A.count(2)
cnt3 = A.count(3)

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













def func(c1, c2, c3):
    if dp[c1][c2][c3] != -1:
        return dp[c1][c2][c3]

    
    prob0 = (N - c1 - c2 - c3) / N

    temp = 1 / (1 - prob0)
    if c1 > 0:
        temp += func(c1 - 1, c2, c3) * (c1 / N) / (1 - prob0)
    if c2 > 0:
        temp += func(c1 + 1, c2 - 1, c3) * (c2 / N) / (1 - prob0)
    if c3 > 0:
        temp += func(c1, c2 + 1, c3 - 1) * (c3 / N) / (1 - prob0)

    dp[c1][c2][c3] = temp
    return temp


print(func(cnt1, cnt2, cnt3))

Solution (Python):



N = int(input())
A = list(map(int, input().split()))
cnt1 = A.count(1)
cnt2 = A.count(2)
cnt3 = A.count(3)

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















def func(c1, c2, c3):
    if dp[c1][c2][c3] != -1:
        return dp[c1][c2][c3]

    
    prob = 1 - (N - c1 - c2 - c3) / N

    temp = 1 / prob
    if c1 > 0:
        temp += func(c1 - 1, c2, c3) * (c1 / N) / prob
    if c2 > 0:
        temp += func(c1 + 1, c2 - 1, c3) * (c2 / N) / prob
    if c3 > 0:
        temp += func(c1, c2 + 1, c3 - 1) * (c3 / N) / prob

    dp[c1][c2][c3] = temp
    return temp


print(func(cnt1, cnt2, cnt3))

Solution (Python):



N = int(input())
A = list(map(int, input().split()))
cnt1 = A.count(1)
cnt2 = A.count(2)
cnt3 = A.count(3)

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













def func(c1, c2, c3):
    if dp[c1][c2][c3] != -1:
        return dp[c1][c2][c3]

    
    prob0 = (N - c1 - c2 - c3) / N

    temp = 1 / (1 - prob0)
    if c1 > 0:
        temp += func(c1 - 1, c2, c3) * (c1 / N) / (1 - prob0)
    if c2 > 0:
        temp += func(c1 + 1, c2 - 1, c3) * (c2 / N) / (1 - prob0)
    if c3 > 0:
        temp += func(c1, c2 + 1, c3 - 1) * (c3 / N) / (1 - prob0)

    dp[c1][c2][c3] = temp
    return temp


print(func(cnt1, cnt2, cnt3))

Solution (Python):
import numpy as np
from numba import njit
 
 
@njit('f8(i8,i8,i8,i8,f8[:,:,:])', cache=True)
def get_expv(N, i, j, k, expv):
    result = 0.0
    if expv[i][j][k] != -1:
        return expv[i][j][k]
    nonlp = 1.0 - (N-i-j-k) / N
    result += 1
    if i != 0:
        result += get_expv(N, i-1, j, k, expv) * (i / N)
    if j != 0:
        result += get_expv(N, i+1, j-1, k, expv) * (j / N)
    if k != 0:
        result += get_expv(N, i, j+1, k-1, expv) * (k / N)
    result /= nonlp
    expv[i][j][k] = result
    return result
 
def main():
    N = int(input())
 
    sushi_count = {1:0, 2:0, 3:0}
    for a in map(int, input().split()):
        sushi_count[a] += 1
    sh = (N+1, N-sushi_count[1]+1, sushi_count[3]+1)
    expv = np.full(sh, -1.0)
 
    expv[0][0][0] = 0.0
 
    print(get_expv(N, sushi_count[1], sushi_count[2], sushi_count[3], expv))
 
if __name__ == '__main__':
    main()

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


@njit('f8(i8,i8,i8,i8,f8[:,:,:])')
def get_expv(N, i, j, k, expv):
    result = 0.0
    if expv[i][j][k] != -1:
        return expv[i][j][k]
    nonlp = 1.0 - (N-i-j-k) / N
    result += 1
    if i != 0:
        result += get_expv(N, i-1, j, k, expv) * (i / N)
    if j != 0:
        result += get_expv(N, i+1, j-1, k, expv) * (j / N)
    if k != 0:
        result += get_expv(N, i, j+1, k-1, expv) * (k / N)
    result /= nonlp
    expv[i][j][k] = result
    return result

def main():
    N = int(input())

    sushi_count = {1:0, 2:0, 3:0}
    for a in map(int, input().split()):
        sushi_count[a] += 1
    sh = (N+1, N-sushi_count[1]+1, sushi_count[3]+1)
    expv = np.full(sh, -1.0)

    expv[0][0][0] = 0.0

    print(get_expv(N, sushi_count[1], sushi_count[2], sushi_count[3], expv))

if __name__ == '__main__':
    main()

Solution (Python):
import sys
input=sys.stdin.readline
def main():
    n=int(input())
    al=tuple(map(int,input().split()))
    arr=[0,0,0]
    for i in al:
        arr[i-1]+=1
    dp=[[[0 for k in range(n+1)] for j in range(n+1)] for i in range(n+1)]
    for c in range(n+1):
        for b in range(n+1-c):
            for a in range(n-b-c+1):
                if a==b==c==0:
                    continue
                if 0<=c-1<=n and 0<=b+1<=n:
                    dp[a][b][c]+=dp[a][b+1][c-1]*c/n
                if 0<=b-1<=n and 0<=a+1<=n:
                    dp[a][b][c]+=dp[a+1][b-1][c]*b/n
                if 0<=a-1<=n:
                    dp[a][b][c]+=dp[a-1][b][c]*a/n
                dp[a][b][c]=(dp[a][b][c]+1)*n/(a+b+c)
    print(dp[arr[0]][arr[1]][arr[2]])
if __name__=='__main__':
    main()

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):
N = int(input())
A = list(map(int, input().split()))


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

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

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

print(mem_rec(A.count(1), A.count(2), A.count(3)))