// 基本テンプレート
 
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <string>
#include <cstring>
#include <deque>
#include <list>
#include <queue>
#include <stack>
#include <vector>
#include <utility>
#include <algorithm>
#include <map>
#include <set>
#include <complex>
#include <cmath>
#include <limits>
#include <cfloat>
#include <climits>
#include <ctime>
#include <cassert>
#include <numeric>
#include <fstream>
#include <functional>
using namespace std;
 
#define rep(i,a,n) for(int (i)=(a); (i)<(n); (i)++)
#define repq(i,a,n) for(int (i)=(a); (i)<=(n); (i)++)
#define repr(i,a,n) for(int (i)=(a); (i)>=(n); (i)--)
#define int long long int
 
template<typename T> void chmax(T &a, T b) {a = max(a, b);}
template<typename T> void chmin(T &a, T b) {a = min(a, b);}
template<typename T> void chadd(T &a, T b) {a = a + b;}
 
typedef pair<int, int> pii;
typedef long long ll;
 
int dx[] = {0, 0, 1, -1};
int dy[] = {1, -1, 0, 0};
const ll INF = 1001001001001001LL;
const ll MOD = 1000000007LL;

int N;
double dp[110][110][110];
int cnt[4];

double solve(int x, int y, int z) {
    if(dp[x][y][z] != -1) return dp[x][y][z];
    dp[x][y][z] = 1;
    if(x != N) dp[x][y][z] += 1.0 * solve(x+1, y, z) * (N - x) / N;
    if(y != x) dp[x][y][z] += 1.0 * solve(x, y+1, z) * (x - y) / N;
    if(z != y) dp[x][y][z] += 1.0 * solve(x, y, z+1) * (y - z) / N;
    dp[x][y][z] *= 1.0 * N / (N - z);
    return dp[x][y][z];
}
 
signed main() {
    cin >> N;
    repq(i,0,N) repq(j,0,N) repq(k,0,N) dp[i][j][k] = -1;
    dp[N][N][N] = 0;

    rep(i,0,N) {
        int A; cin >> A;
        cnt[min(A, 3LL)]++;
    }
    repr(i,2,0) {
        cnt[i] += cnt[i+1];
    }

    printf("%.12f\n", solve(cnt[1], cnt[2], cnt[3]));
    return 0;
}