Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
sys.setrecursionlimit(10 ** 8)
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1

divT = [1/i for i in range(1,n+1)]

from functools import lru_cache### 
@lru_cache(100100100)

def f(x,y,z):
    if x==y==z==n:
        return 0
    res = 0
    delta = n*divT[n-z-1]
    if x != n:
        res += (f(x+1,y,z) + delta) * (n-x)
    if y != x:
        res += (f(x,y+1,z) + delta) * (x-y)#x+1,y+1ではない
    if z != y:
        res += (f(x,y,z+1) + delta)* (y-z)#x+1,y+1,z+1ではない
    res *=divT[n-z-1]
    return res
    

print(f(cnt[1]+cnt[2],cnt[2],0))


Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
sys.setrecursionlimit(10 ** 8)
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1

divT = [1/i for i in range(1,n+1)]

from functools import lru_cache### 
@lru_cache(100100100)

def f(x,y,z):
    if x==y==z==n:
        return 0
    res = 0
    delta = n*divT[n-z-1]
    if x != n:
        res += (f(x+1,y,z) + delta) * (n-x)
    if y != x:
        res += (f(x,y+1,z) + delta) * (x-y)#x+1,y+1ではない
    if z != y:
        res += (f(x,y,z+1) + delta)* (y-z)#x+1,y+1,z+1ではない
    res *=divT[n-z-1]
    return res
    

print(f(cnt[1]+cnt[2],cnt[2],0))


Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Runtime Error):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1

divT = [1/i for i in range(1,n+1)]

from functools import lru_cache### 
@lru_cache(100100100)

def f(x,y,z):
    if x==y==z==n:
        return 0
    res = 0
    delta = n*divT[n-z-1]
    if x != n:
        res += (f(x+1,y,z) + delta) * (n-x)
    if y != x:
        res += (f(x,y+1,z) + delta) * (x-y)#x+1,y+1ではない
    if z != y:
        res += (f(x,y,z+1) + delta)* (y-z)#x+1,y+1,z+1ではない
    res *=divT[n-z-1]
    return res
    

print(f(cnt[1]+cnt[2],cnt[2],0))


Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Runtime Error):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1

divT = [1/i for i in range(1,n+1)]

from functools import lru_cache### これの練習
@lru_cache(100000)

def f(x,y,z):
    if x==y==z==n:
        return 0
    res = 0
    delta = n*divT[n-z-1]
    if x != n:
        res += (f(x+1,y,z) + delta) * (n-x)
    if y != x:
        res += (f(x,y+1,z) + delta) * (x-y)#x+1,y+1ではない
    if z != y:
        res += (f(x,y,z+1) + delta)* (y-z)#x+1,y+1,z+1ではない
    res *=divT[n-z-1]
    return res
    

print(f(cnt[1]+cnt[2],cnt[2],0))


Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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


Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[-1 for _ in range(n+2)] for __ in range(n+2)] for ___ in range(n+2)]
"""

for i in range(n,-1,-1):
    for j in range(n,i-1,-1):
        for k in range(n,j-1,-1):
            dp[k][j][i] = (dp[k-1][j+1][i] + 1) * k/n +(dp[k][j-1][i+1] + 1) * i/n +  (dp[k][j][i] + 1) * j/n
こういうインデックスが面倒なことになるから1枚、2枚、3枚でなく1枚以上、2枚以上、3枚以上食べた種類の数で管理する
"""
"""
for i in range(n-1,-1,-1):
    for j in range(n-1,i-1,-1):
        for k in range(n-1,j-1,-1):
            
            ???
            dp[k][j][i] = (dp[k+1][j+1][i] + 1) * k/n +(dp[k+1][j+1][i+1] + 1) * j/n +  (dp[k][j][i] + 1) * i/n
            → dp[k][j][i] = ((dp[k+1][j+1][i] + 1) * k/(n-i) +(dp[k+1][j+1][i+1] + 1) * j/(n-i)) +   n/i
            """


#上の方針だと追えないのでやり直し。
#割り算が遅いと聞いてとりあえず試す
dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]
def solve(x, y, z):
    if dp[x][y][z] != -1:
        return dp[x][y][z]
    else:
        delta = n*divT[n-z-1]#3枚以上ある皿以外を引くまで
        dp[x][y][z] = 0
        if x != n:
            dp[x][y][z] += (solve(x+1,y,z) + delta) * (n-x)*divT[n-z-1]
        if y != x:
            dp[x][y][z] += (solve(x,y+1,z) + delta) * (x-y)*divT[n-z-1]#x+1,y+1ではない
        if z != y:
            dp[x][y][z] += (solve(x,y,z+1) + delta) * (y-z)*divT[n-z-1]#x+1,y+1,z+1ではない
       
        return dp[x][y][z]


print(solve(cnt[1]+cnt[2], cnt[2] ,0))




Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[-1 for _ in range(n+2)] for __ in range(n+2)] for ___ in range(n+2)]
"""

for i in range(n,-1,-1):
    for j in range(n,i-1,-1):
        for k in range(n,j-1,-1):
            dp[k][j][i] = (dp[k-1][j+1][i] + 1) * k/n +(dp[k][j-1][i+1] + 1) * i/n +  (dp[k][j][i] + 1) * j/n
こういうインデックスが面倒なことになるから1枚、2枚、3枚でなく1枚以上、2枚以上、3枚以上食べた種類の数で管理する
"""
"""
for i in range(n-1,-1,-1):
    for j in range(n-1,i-1,-1):
        for k in range(n-1,j-1,-1):
            
            ???
            dp[k][j][i] = (dp[k+1][j+1][i] + 1) * k/n +(dp[k+1][j+1][i+1] + 1) * j/n +  (dp[k][j][i] + 1) * i/n
            → dp[k][j][i] = ((dp[k+1][j+1][i] + 1) * k/(n-i) +(dp[k+1][j+1][i+1] + 1) * j/(n-i)) +   n/i
            """


#上の方針だと追えないのでやり直し。
#割り算が遅いと聞いてとりあえず試す
dp[n][n][n] = 0
def solve(x, y, z):
    if dp[x][y][z] != -1:
        return dp[x][y][z]
    else:
        
        hoge = 1.0 / (n - z)
        delta = n*hoge#3枚以上ある皿以外を引くまで
        dp[x][y][z] = 0
        try:
            if x != n:
                dp[x][y][z] += (solve(x+1,y,z) + delta) * (n-x)*hoge
            if y != x:
                dp[x][y][z] += (solve(x,y+1,z) + delta) * (x-y)*hoge#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (solve(x,y,z+1) + delta) * (y-z)*hoge#x+1,y+1,z+1ではない
        except:
            print(x,y,z)

        return dp[x][y][z]


print(solve(cnt[1]+cnt[2], cnt[2] ,0))




Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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

Buggy python submission (verdict: Time Limit Exceeded):
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
n = int(input())
a = list(map(int, input().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[-1 for _ in range(n+2)] for __ in range(n+2)] for ___ in range(n+2)]
"""

for i in range(n,-1,-1):
    for j in range(n,i-1,-1):
        for k in range(n,j-1,-1):
            dp[k][j][i] = (dp[k-1][j+1][i] + 1) * k/n +(dp[k][j-1][i+1] + 1) * i/n +  (dp[k][j][i] + 1) * j/n
こういうインデックスが面倒なことになるから1枚、2枚、3枚でなく1枚以上、2枚以上、3枚以上食べた種類の数で管理する
"""
"""
for i in range(n-1,-1,-1):
    for j in range(n-1,i-1,-1):
        for k in range(n-1,j-1,-1):
            
            ???
            dp[k][j][i] = (dp[k+1][j+1][i] + 1) * k/n +(dp[k+1][j+1][i+1] + 1) * j/n +  (dp[k][j][i] + 1) * i/n
            → dp[k][j][i] = ((dp[k+1][j+1][i] + 1) * k/(n-i) +(dp[k+1][j+1][i+1] + 1) * j/(n-i)) +   n/i
            """


#上の方針だと追えないのでやり直し。
dp[n][n][n] = 0
def solve(x, y, z):
    if dp[x][y][z] != -1:
        return dp[x][y][z]
    else:
        delta = n/(n-z)#3枚以上ある皿以外を引くまで
        dp[x][y][z] = 0
        if x != n:
            dp[x][y][z] += (solve(x+1,y,z) + delta) * (n-x)/(n-z)
        if y != x:
                dp[x][y][z] += (solve(x,y+1,z) + delta) * (x-y)/(n-z)#x+1,y+1ではない
        if z != y:
            dp[x][y][z] += (solve(x,y,z+1) + delta) * (y-z)/(n-z)#x+1,y+1,z+1ではない

        return dp[x][y][z]


print(solve(cnt[1]+cnt[2], cnt[2] ,0))




Fixed submission:
#http://compro.tsutajiro.com/archive/180220_probability_dp.pdf
import sys
n = int(input())
a = list(map(int, sys.stdin.readline().split( )))
cnt =[0,0,0]
for i in range(n):
    cnt[3-a[i]] += 1
dp = [[[0 for _ in range(n+1)] for __ in range(n+1)] for ___ in range(n+1)]

dp[n][n][n] = 0
divT = [1/i for i in range(1,n+1)]

for x in range(n,cnt[1]+cnt[2]-1,-1):## -2だったのを-1に
    for y in range(x,cnt[2]-1,-1):##cnt[2]
        for z in range(y,-1,-1):
            delta = n*divT[n-z-1]#2枚以下を引くまで
            if x != n:
                dp[x][y][z] += (dp[x+1][y][z] + delta) * (n-x)
            if y != x:
                dp[x][y][z] += (dp[x][y+1][z] + delta) * (x-y)#x+1,y+1ではない
            if z != y:
                dp[x][y][z] += (dp[x][y][z+1] + delta)* (y-z)#x+1,y+1,z+1ではない
            dp[x][y][z] *=divT[n-z-1]

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