Buggy python submission (verdict: Time Limit Exceeded):
#事象が確率pで起きるとき, 事象が起きるまでの回数の期待値は1/p
#https://qiita.com/drken/items/03c7db44ccd27820ea0d

def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10000000)
    from collections import Counter, deque
    #from collections import defaultdict
    from itertools import combinations, permutations, accumulate
    #from itertools import product
    from bisect import bisect_left,bisect_right
    import heapq
    from math import floor, ceil
    #from operator import itemgetter

    #inf = 10**17
    #mod = 10**9 + 7

    N = int(input())
    a = tuple(map(int, input().split()))
    c1 = a.count(1) #i
    c2 = a.count(2) #j
    c3 = a.count(3) #k

    dp = [[[0]*(N+1) for _ in range(N+1)] for __ in range(N+1)]
    for i in range(c3+1):
        for j in range(c2+c3-i+1):
            for k in range(N-i-j+1):
                if i == j == k == 0:
                    continue
                tmp = N
                if k != 0:
                    tmp += k*dp[k-1][j][i]
                if j != 0:
                    tmp += j*dp[k+1][j-1][i]
                if i != 0:
                    tmp += i*dp[k][j+1][i-1]
                dp[k][j][i] = tmp/(i+j+k)

    print(dp[c1][c2][c3])

if __name__ == '__main__':
    main()

Fixed submission:
def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10000000)
    from collections import Counter, deque
    #from collections import defaultdict
    from itertools import combinations, permutations, accumulate
    #from itertools import product
    from bisect import bisect_left,bisect_right
    import heapq
    from math import floor, ceil
    #from operator import itemgetter

    #inf = 10**17
    #mod = 10**9 + 7

    N = int(input())
    a = tuple(map(int, input().split()))
    c1 = a.count(1) #i
    c2 = a.count(2) #j
    c3 = a.count(3) #k

    dp = [[[0.0]*(N+1) for _ in range(N+1)] for __ in range(N+1)]
    for i in range(c3+1):
        for j in range(c2+c3-i+1):
            for k in range(N-i-j+1):
                if i == j == k == 0:
                    continue
                tmp = N
                if k != 0:
                    tmp += k*dp[k-1][j][i]
                if j != 0:
                    tmp += j*dp[k+1][j-1][i]
                if i != 0:
                    tmp += i*dp[k][j+1][i-1]
                dp[k][j][i] = tmp/(i+j+k)

    print(dp[c1][c2][c3])

if __name__ == '__main__':
    main()

Buggy python submission (verdict: Time Limit Exceeded):
#事象が確率pで起きるとき, 事象が起きるまでの回数の期待値は1/p
#https://qiita.com/drken/items/03c7db44ccd27820ea0d

def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10000000)
    from collections import Counter, deque
    #from collections import defaultdict
    from itertools import combinations, permutations, accumulate
    #from itertools import product
    from bisect import bisect_left,bisect_right
    import heapq
    from math import floor, ceil
    #from operator import itemgetter

    #inf = 10**17
    #mod = 10**9 + 7

    N = int(input())
    a = tuple(map(int, input().split()))
    c1 = a.count(1) #i
    c2 = a.count(2) #j
    c3 = a.count(3) #k

    dp = [[[0]*(N+1) for _ in range(N+1)] for __ in range(N+1)]
    for i in range(c3+1):
        for j in range(c2+c3-i+1):
            for k in range(N-i-j+1):
                if i == j == k == 0:
                    continue
                tmp = N
                if k != 0:
                    tmp += k*dp[k-1][j][i]
                if j != 0:
                    tmp += j*dp[k+1][j-1][i]
                if i != 0:
                    tmp += i*dp[k][j+1][i-1]
                dp[k][j][i] = tmp/(i+j+k)

    print(dp[c1][c2][c3])

if __name__ == '__main__':
    main()

Fixed submission:
def main():
    import sys
    input = sys.stdin.readline
    sys.setrecursionlimit(10000000)
    from collections import Counter, deque
    #from collections import defaultdict
    from itertools import combinations, permutations, accumulate
    #from itertools import product
    from bisect import bisect_left,bisect_right
    import heapq
    from math import floor, ceil
    #from operator import itemgetter

    #inf = 10**17
    #mod = 10**9 + 7

    N = int(input())
    a = tuple(map(int, input().split()))
    c1 = a.count(1) #i
    c2 = a.count(2) #j
    c3 = a.count(3) #k

    dp = [[[0.0]*(N+1) for _ in range(N+1)] for __ in range(N+1)]
    for i in range(c3+1):
        for j in range(c2+c3-i+1):
            for k in range(N-i-j+1):
                if i == j == k == 0:
                    continue
                tmp = N
                if k != 0:
                    tmp += k*dp[k-1][j][i]
                if j != 0:
                    tmp += j*dp[k+1][j-1][i]
                if i != 0:
                    tmp += i*dp[k][j+1][i-1]
                dp[k][j][i] = tmp/(i+j+k)

    print(dp[c1][c2][c3])

if __name__ == '__main__':
    main()