#include <bits/stdc++.h>
#define elif else if
#define ALL(a) a.begin(), a.end()
using namespace std;
using LL = long long;
using LD = long double;
using PII = pair<int, int>;
using PLI = pair<long long, int>;
using PLL = pair<long long, long long>;
using TI = tuple<int, int, int>;
mt19937 mrand(random_device{}());
int rnd(int x) { return mrand() % x; }
template <class T> int bcnt(const T &a) { return __builtin_popcount(a); }
template <> int bcnt(const int64_t &a) { return __builtin_popcountll(a); }
template <class T> int blen(const T &a) { return a ? 1 + blen(a / 2) : 0; }
template <class A, class B> B conv(const A &a, B b = B()) {
    stringstream s;
    s << a, s >> b;
    return b;
}
template <class T> T sign(const T &a) { return a == 0 ? 0 : (a < 0 ? -1 : 1); }
template <class T> T floor(const T &a, const T &b) {
    assert(b != 0);
    return sign(a) * sign(b) > 0 ? abs(a) / abs(b) : -(abs(a) + abs(b) - 1) / abs(b);
}
template <class T> T ceil(const T &a, const T &b) {
    assert(b != 0);
    return sign(a) * sign(b) > 0 ? (abs(a) + abs(b) - 1) / abs(b) : -abs(a) / abs(b);
}
template <class T> bool chmin(T &a, T b) { return b < a ? a = b, true : false; }
template <class T> bool chmax(T &a, T b) { return b > a ? a = b, true : false; }
template <class T> auto min(const T &a) { return *min_element(ALL(a)); } 
template <class T> auto max(const T &a) { return *max_element(ALL(a)); }
template <class T> auto sum(const T &a) { return accumulate(ALL(a), (typename T::value_type)0); }
template <> auto sum(const vector<int> &a) { return accumulate(ALL(a), 0LL); }
template <> auto sum(const vector<string> &a) { return accumulate(ALL(a), string()); }
template <class T> void uniq(T &a) { sort(ALL(a)); a.erase(unique(ALL(a)), a.end()); }
template<class A, class B> ostream& operator <<(ostream& out, const pair<A, B> &p) {
    return out << "(" << p.first << ", " << p.second << ")";
}
template <typename T, size_t N>
ostream& operator <<(ostream& out, const array<T, N>& a) {
    out << "["; bool first = true;
    for (auto& v : a) { out << (first ? "" : ", "); out << v; first = 0;} out << "]";
    return out;
}
template <typename T>
ostream& operator <<(ostream& out, const vector<T>& a) {
    out << "["; bool first = true;
    for (auto& v : a) { out << (first ? "" : ", "); out << v; first = 0;} out << "]";
    return out;
}
template <typename T>
ostream& operator <<(ostream& out, const deque<T>& a) {
    out << "["; bool first = true;
    for (auto& v : a) { out << (first ? "" : ", "); out << v; first = 0;} out << "]";
    return out;
}
template <typename T, class Cmp>
ostream& operator <<(ostream& out, const set<T, Cmp>& a) {
    out << "{"; bool first = true;
    for (auto& v : a) { out << (first ? "" : ", "); out << v; first = 0;} out << "}";
    return out;
}
template <typename T>
ostream& operator <<(ostream& out, const unordered_set<T>& a) {
    out << "{"; bool first = true;
    for (auto& v : a) { out << (first ? "" : ", "); out << v; first = 0;} out << "}";
    return out;
}
template <typename U, typename T, class Cmp>
ostream& operator <<(ostream& out, const map<U, T, Cmp>& a) {
    out << "{"; bool first = true;
    for (auto& p : a) { out << (first ? "" : ", "); out << p.first << ":" << p.second; first = 0;} out << "}";
    return out;
}
template<class T> inline void YES(T condition) { 
    if(condition) cout << "YES" << '\n';
    else cout << "NO" << '\n';
}
template<class T> inline void Yes(T condition) {
    if(condition) cout << "Yes" << '\n';
    else cout << "No" << '\n';
}
template<class T> T gcd(T a, T b) {
    while (b) {
        T tmp = b;
        b = a % b;
        a = tmp;
    }
    return a;
}
template<class T> T lcm(T a, T b){
    return a / gcd(a, b) * b;
}

#define trace(...) __f(#__VA_ARGS__, __VA_ARGS__)
template <typename Arg1>
void __f(const char* name, Arg1&& arg1){
    cerr << name << ": " << arg1 << endl;
}
template <typename Arg1, typename... Args>
void __f(const char* names, Arg1&& arg1, Args&&... args){
    const char* comma = strchr(names + 1, ',');
    cerr.write(names, comma - names) << ": " << arg1 << " |";
    __f(comma + 1, args...);
}

int qpow(int a, int b, int MOD) {
    return qpow(static_cast<long long>(a), static_cast<long long>(b), MOD);
}
long long qpow(long long a, long long b, int MOD) {
    assert(b >= 0);
    if (b == 0) return 1;
    long long res = 1;
    while (b) {
        if (b & 1) res = (res * a) % MOD;
        a = (a * a) % MOD;
        b >>= 1;
    }
    return res;
}
template <class T> auto vect(const T &v, int n) { return vector<T>(n, v); }
template <class T, class... D> auto vect(const T& v, int n, D... m) {
    return vector<decltype(vect(v, m...))>(n, vect(v, m...));
}

// template <typename T>
// vector<vector<T> > matmul(vector<vector<T> >& a, vector<vector<T> >& b) {
//     int n = a.size(), m = b.size();
//     assert(n > 0 && m > 0);
//     assert((int)a[0].size() == m);
//     int p = b[0].size();
//     vector<vector<T> > res(n, vector<T>(p, (typename T::value_type)0));
//     for (int i = 0; i < n; i++) {
//         for (int k = 0; k < p; k++) {
//             T now = (typename T::value_type)0;
//             for (int j = 0; j < m; j++) {
//                 now += a[i][j] * b[j][k];
//             }
//             res[i][k] = now;
//         }
//     }
// }
// template <>
// vector<vector<LL> > matmul(vector<vector<LL> >& a, vector<vector<LL> >& b, int MOD) {
//     int n = a.size(), m = b.size();
//     assert(n > 0 && m > 0);
//     assert((int)a[0].size() == m);
//     int p = b[0].size();
//     vector<vector<LL> > res(n, vector<LL>(p, 0LL));
//     for (int i = 0; i < n; i++) {
//         for (int k = 0; k < p; k++) {
//             LL now = (typename T::value_type)0;
//             for (int j = 0; j < m; j++) {
//                 now += a[i][j] * b[j][k];
//             }
//             res[i][k] = now % MOD;
//         }
//     }
// }
// vector<vector<LL> > qpow(vector<vector<LL> >& a, long long b, int MOD) {
//     assert(b >= 0);
//     int n = a.size();
//     assert(n >= 1);
//     int m = a[0].size();
//     assert(n == m);
//     if (b == 0) {
//         vector<vector<LL> > res(n, vector<LL>(n, 0LL));
//         for (int i = 0; i < n; i++) res[i][i] = 1LL;
//         return res;
//     }
//     while (b) {
//         if (b & 1) 
//     }

// }


// long long qpow(long long a, long long b, int MOD) {
//     if (b == 0) return 1;
//     long long res = 1;
//     while (b) {
//         if (b & 1) res = (res * a) % MOD;
//         a = (a * a) % MOD;
//         b >>= 1;
//     }
//     return res;
// }

const int MOD = 1e9 + 7;
const int iinf = __INT_MAX__ / 2;
const LL linf = __INT64_MAX__ / 2;
const int maxn = 2e6;
const LD eps = 1e-8;
const LD pi = acos(-1);
// const int INF = 2e9;
int fac[maxn + 5];
int invfac[maxn + 5];
void init() {
    fac[0] = 1;
    for (int i = 1; i <= maxn; i++) fac[i] = 1LL * fac[i - 1] * i % MOD;
    invfac[maxn] = (int)(qpow(fac[maxn], MOD - 2, MOD));
    for (int i = maxn - 1; i >= 0; i--) invfac[i] = 1LL * invfac[i + 1] * (i + 1) % MOD;
}
int C(int n, int k) {
    assert(0 <= k && n >= k);
    int tmp = (int)(1LL * invfac[k] * invfac[n - k] % MOD);
    return (int)(1LL * fac[n] * tmp % MOD);
}
    // vector<int> fa(n);
    // vector<int> sz(n, 1);
    // iota(fa.begin(), fa.end(), 0);
    // function<int(int)> Find = [&](int x) {
    //     return x == fa[x] ? x : fa[x] = Find(fa[x]);
    // };
    // function<bool(int, int)> Union = [&](int x, int y) {
    //     int px = Find(x), py = Find(y);
    //     if (px == py) return false;
    //     if (sz[px] > sz[py]) swap(px, py);
    //     fa[px] = py;
    //     sz[py] += sz[px];
    //     return true;
    // };

int main() {
    cin.tie(nullptr);
    ios::sync_with_stdio(false);
    cout << fixed << setprecision(20);

    int n; cin >> n;
    vector<int> cnts(4);
    for (int i = 0; i < n; i++) {
        int x; cin >> x;
        cnts[x]++;
    }
    auto dp = vect<double>(-1.0, n + 1, n + 1, n + 1);
    function<double(int, int, int)> dfs = [&](int c1, int c2, int c3) {
        if (!c1 && !c2 && !c3) return 0.0;
        double& res = dp[c1][c2][c3];
        if (res >= 0.0) return res;
        int tot = c1 + c2 + c3;
        int waste = n - c1 - c2 - c3;
        double wasteP = static_cast<double>(waste) / static_cast<double>(n);
        // expectation value of waste
        // p = usefulP = tot / n
        // Evw = p(1 - p) + 2p(1 - p)^2 + 3p(1-p)^3 + ...
        double Evw = static_cast<double>(n - tot) / static_cast<double>(tot);
        double Ev1 = 0.0, Ev2 = 0.0, Ev3 = 0.0;
        if (c1) Ev1 = (dfs(c1 - 1, c2, c3) + 1.0) * c1 / tot;
        if (c2) Ev2 = (dfs(c1 + 1, c2 - 1, c3) + 1.0) * c2 / tot;
        if (c3) Ev3 = (dfs(c1, c2 + 1, c3 - 1) + 1.0) * c3 / tot;
        // trace(c1, c2, c3, Evw, Ev1, Ev2, Ev3);
        return res = Evw + Ev1 + Ev2 + Ev3;
    };
    cout << dfs(cnts[1], cnts[2], cnts[3]) << '\n';

    return 0;
}


Measured memory_bytes: 219920