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



Sample input:
3
1 1 1


Sample output:
5.5

Submission status: Accepted

Solution code:
# 期待値自体が苦手なので厳しい
# 解説AC
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()

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
# ===============================================================================================
# importing some useful libraries.
from __future__ import division, print_function
from fractions import Fraction
import sys
import os
from io import BytesIO, IOBase

# If the element is already present in the list,

# the left most position where element has to be inserted is returned.
from bisect import bisect_right as br
from bisect import bisect

# If the element is already present in the list,
# the right most position where element has to be inserted is returned

# ==============================================================================================
# fast I/O region

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)

# inp = lambda: sys.stdin.readline().rstrip("\r\n")


# ===============================================================================================
# some shortcuts

mod = 1000000007


def inp(): return sys.stdin.readline().rstrip("\r\n")  # for fast input


def out(var): sys.stdout.write(str(var))  # for fast output, always take string


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")  # as stdout.write always print sring.


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]])

Submission status: Accepted

Solution code:
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))
# 受け渡されたすべての要素からsubだけ引く.リストを*をつけて展開しておくこと


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)  # byte同士の結合ができないのでreadlineでなくreadで
    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  # 2147483648 > 10**9
# default import

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  # , insort_left, insort_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[:])',cache=True)


@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))

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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))
# 受け渡されたすべての要素からsubだけ引く.リストを*をつけて展開しておくこと


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)  # byte同士の結合ができないのでreadlineでなくreadで
    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  # 2147483648 > 10**9
# default import

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  # , insort_left, insort_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[:])',cache=True)


@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))

Submission status: Accepted

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

# sys.setrecursionlimit(10**4)
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]])

Submission status: Accepted

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

Submission status: Accepted

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

# PyPyでは事前の型推論によりチェックを一部省略することで高速化を図っている部分があるが、
# 処理の中でintもfloatも取り得る変数にしてしまうとその推論が上手く働かない。
# 今回のように小数が入るとわかっている場合は、初期化する値も小数にしておけば、よりPyPyの恩恵を享受できる
# ↑ほんと？

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))

Submission status: Accepted

Solution code:
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))

Submission status: Accepted

Solution code:
# dp[c1][c2][c3]=
# 		(1個の皿がc1枚、2個の皿がc2枚、3個の皿がc3枚あるとき、全て食べるまでの操作回数の期待値)

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

# dp[c1][c2][c3]=1
# 	+dp[c1−1][c2][c3]∗(1個の皿が選ばれる確率)
# 	+dp[c1+1][c2−1][c3]∗(2個の皿が選ばれる確率)
# 	+dp[c1][c2+1][c3−1]∗(3個の皿が選ばれる確率)
# 	+dp[c1][c2][c3]∗(0個の皿が選ばれる確率)
# ↓式変形↓
# dp[c1][c2][c3]=1/(1−(0個の皿が選ばれる確率))
# 			+ dp[c1−1][c2][c3]*(1個の皿が選ばれる確率)/(1−(0個の皿が選ばれる確率))
# 			+ dp[c1+1][c2−1][c3]*(2個の皿が選ばれる確率)/(1−(0個の皿が選ばれる確率))
# 			+ dp[c1][c2+1][c3−1]*(3個の皿が選ばれる確率)/(1−(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))

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
import sys
sys.setrecursionlimit(10**7)
input = sys.stdin.readline
def main():
    n=int(input())
    a=list(map(int,input().split()))
    # 0個の皿がx0枚、1個の皿がx1枚、2個の皿がx2枚、3個の皿がx3枚
    # x0=n-x1-x2-x3
    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()

Submission status: Accepted

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

# 残り１、２、３枚の皿がi, j, k種類ある時に食べきるのに必要な回数をdpで求める
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)))

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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))

Submission status: Accepted

Solution code:
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)