\(\mathcal{Description}\)

  Link.

  给定非负整数序列 \(\{a_n\}\) 和 \(m\)，每次随机在 \(\{a\}\) 中取一个非零的 \(a_i\)（保证存在），令其 \(-1\),重复 \(m\) 次，求最终 \(\{a\}\) 中 \(0\) 的期望个数。

  \(n\le15\)，\(m\le100\)。

\(\mathcal{Solution}\)

  既然有 DP of DP，我觉得这种题就叫 DP and DP。（

  看到 \(n\le15\)，想要状压，但无论如何都不能完整记录 \(\{a\}\) 的状态，也就不可能及时知道某个 \(a\) 会否变成 \(0\)，怎么办呢？

  既然不知道，我们索性不去动它们——设前 \(i\) 次操作后，\(\{a\}\) 中为 \(0\) 的下标集合为 \(S\)， 我们只落实对 \(S\) 中元素造成影响的操作，共 \(\sum_{i\in S}a_i\) 次；而对于剩下 \(\left(i-\sum_{i\in S}a_i\right)\) 次，只记录为“它们会影响某个不在 \(T\) 中元素”，其中 \(T\) 是进行这个操作时的 \(0\) 值集合；“某个”是具体的，但尚未确定的一个，所以提供概率 \(\frac{1}{n-|T|}\)。另一方面，只有在扩充 \(S\) 时来确定若干“某个”，此时就只需要组合选取，而不需要乘概率了。那么摆出式子，令 \(f(i,S)\) 表示 \(i\) 次操作后，\(0\) 值集合为 \(S\) 的概率（概率计算方式如前文），则：

\[f(i,S)=\frac{1}{n-|S|}f(i-1,S)+\frac{1}{n-|S|+1}\sum_{j\in S}\binom{i-1-\sum_{k\in S,k\not= j}a_k}{a_j-1}f(i-1,S\setminus\{j\}).
\]

  思考此时得到的，应该用于计算答案的 \(f(m,S)\) 意味什么：落在 \(S\) 内的操作，它们是正常的，而其余的 \(\left(m-\sum_{i\in S}a_i\right)\) 次，它们一定落在 \(U\setminus S\)（\(U\) 为全集）中，遗漏的仅仅是方案数！所以在 \(f\) 的目标状态为基础，再来一个 DP：令 \(g(i,S)\) 表示有 \(i\) 次操作落在 \(S\)，且 \(S\) 中不存在 \(0\) 值的方案数，简单地，任取一个 \(x\in S\)，得到转移：

\[g(i,S)=\sum_{j=0}^{\min\{i,a_{x}-1\}}\binom{i}{j}g(i-j,S\setminus\{x\}).
\]

最终答案则为：

\[\sum_{S\subseteq U}f(m,S)\cdot g(m-\sum_{i\in S}a_i,U\setminus S)\cdot|S|.
\]

  复杂度为 \(\mathcal O(m(n+m)2^n)\)。

/*~Rainybunny~*/

#include <cstdio>

#define rep( i, l, r ) for ( int i = l, rep##i = r; i <= rep##i; ++i )
#define per( i, r, l ) for ( int i = r, per##i = l; i >= per##i; --i )

inline int imin( const int a, const int b ) { return a < b ? a : b; }

typedef long double LD;
#define double LD

const int MAXN = 15, MAXM = 100;
int n, m, a[MAXN + 5], sum[1 << MAXN], bitc[1 << MAXN], bitw[1 << MAXN];
double comb[MAXM + 5][MAXM + 5], f[2][1 << MAXN], g[MAXM + 5][1 << MAXN];

inline void init() {
    comb[0][0] = 1.;
    rep ( i, 1, m ) {
        comb[i][0] = 1.;
        rep ( j, 1, i ) comb[i][j] = comb[i - 1][j] + comb[i - 1][j - 1];
    }
}

inline void getF() {
    f[0][0] = 1.;
    rep ( i, 1, ( 1 << n ) - 1 ) bitc[i] = bitc[i ^ ( i & -i )] + 1;
    for ( int i = 1, sta = 1; i <= m; ++i, sta ^= 1 ) {
        rep ( S, 0, ( 1 << n ) - 1 ) {
            f[sta][S] = 0.;
            rep ( j, 0, n - 1 ) if ( S >> j & 1 ) {
                f[sta][S] += i < sum[S ^ 1 << j] + 1 ? 0 :
                  comb[i - sum[S ^ 1 << j] - 1][a[j] - 1] * f[!sta][S ^ 1 <<j];
            }
            f[sta][S] /= n - bitc[S] + 1;
            if ( n > bitc[S] ) f[sta][S] += f[!sta][S] / ( n - bitc[S] );
            // fprintf( stderr, "f(%d,%d)=%f\n", i, S, f[sta][S] );
        }
    }
}

inline void getG() {
    rep ( i, 2, ( 1 << n ) - 1 ) bitw[i] = bitw[i >> 1] + 1;
    rep ( S, 0, ( 1 << n ) - 1 ) g[0][S] = 1.;
    rep ( i, 1, m ) {
        rep ( S, 1, ( 1 << n ) - 1 ) {
            int v = bitw[S & -S]; double &cur = g[i][S];
            rep ( j, 0, imin( i, a[v] - 1 ) ) {
                cur += g[i - j][S ^ 1 << v] * comb[i][j];
            }
            // fprintf( stderr, "g(%d,%d)=%.0f\n", i, S, g[i][S] );
        }
    }
}

int main() {
    scanf( "%d %d", &n, &m ), init();
    rep ( i, 0, n - 1 ) scanf( "%d", &a[i] ), sum[1 << i] = a[i];
    rep ( S, 1, ( 1 << n ) - 1 ) sum[S] = sum[S & -S] + sum[S ^ ( S & -S )];

    getF(), getG();

    double ans = 0.;
    rep ( S, 0, ( 1 << n ) - 1 ) if ( m >= sum[S] ) {
        ans += f[m & 1][S] * g[m - sum[S]][repS ^ S] * bitc[S];
    }
    printf( "%.12Lf\n", ans );
    return 0;
}