「集训队作业2018」喂鸽子

设 \(F(n)\) 表示有 \(n\) 只鸽子，每次等概率选一只喂，期望喂饱第一只鸽子的时间，\(f_{n,m}\) 表示有 \(n\) 只鸽子，已经喂了 \(m\) 次，此时这 \(n\) 只鸽子中没有鸽子被喂饱的概率。
\[ Ans = \sum_{i=1}^n (-1)^{i+1}{n \choose i} F(i) \\ F(n)=\sum_{i\geq 0}\sum_{j \leq i}{i\choose j}f_{n,j} (\frac{N-n}{N})^{i-j}\\ =\sum_{i\geq 0}f_{n,i}\sum_{j \geq 0} {i+j\choose i}(\frac{N-n}{N})^{j} \]

注意到有
\[ (\dfrac{1}{1-x})^n=\sum_{i\geq 0} {n+i-1\choose n-1}x^i \]
所以
\[ F(n)=\sum_{i\geq 0}f_{n,i}\sum_{j \geq 0} {i+j\choose i}(\frac{N-n}{N})^{j} \\ =\sum_{i \geq 0}f_{n,i}(\frac{N}{n})^{i+1} \]
题外话：我们要求的实际上是恰有 \(i\) 次分配到这个集合的所有方案里分配到集合外的概率之和，不太好组合意义理解，如果有大爷会组合意义了教教我啊。

接下来求 \(f_{n,i}\) ，这里由于之前没有考虑分配 \(n\) 中玉米的概率，所以在这里最后一只鸽子分配到一个玉米的概率是 \(\frac{1}{N}\) 。
\[ f_{n,i}=\frac{1}{N^i}[x^i] (\sum_{j=0}^{k-1}\frac{x^i}{i!})^n \]
或者也可以递推
\[ f_{n,i} =\sum_{j=0}^{i}{i\choose j}f_{n-1,i-j}\frac{1}{N^j} \]
这样做复杂度是 \(\mathcal O(n^2klog(nk))\) 的，还有复杂度为 \(\mathcal O(n^2k)\) 的高妙做法，以后再补。

code

/*program by mangoyang*/
#include <bits/stdc++.h>
#define inf (0x7f7f7f7f)
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
    int ch = 0, f = 0; x = 0;
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
    for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    if(f) x = -x;
}
const int N = 55, K = 1005, mod = 998244353, G = 3;
int f[N][1<<17], js[1<<17], a[1<<17], inv[1<<17], n, k, ans;
inline void up(int &x, int y){
    x = x + y >= mod ? x + y - mod : x + y;
}
inline int Pow(int a, int b){
    int ans = 1;
    for(; b; b >>= 1, a = 1ll * a * a % mod)
        if(b & 1) ans = 1ll * ans * a % mod;
    return ans;
}
namespace poly{
    int rev[1<<17], len, lg;
    inline void timesinit(int lenth){
        for(len = 1, lg = 0; len <= lenth; len <<= 1, lg++);
        for(int i = 0; i < len; i++)
            rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (lg - 1)); 
    }
    inline void dft(int *a, int sgn){
        for(int i = 0; i < len; i++)
            if(i < rev[i]) swap(a[i], a[rev[i]]);
        for(int k = 2; k <= len; k <<= 1){
            int w = Pow(G, (mod - 1) / k);
            if(sgn == -1) w = Pow(w, mod - 2);
            for(int i = 0; i < len; i += k){
                int now = 1;
                for(int j  = i; j < i + (k >> 1); j++){
                    int x = a[j], y = 1ll * a[j+(k>>1)] * now % mod;
                    a[j] = x + y >= mod ? x + y - mod : x + y;
                    a[j+(k>>1)] = x - y < 0 ? x - y + mod : x - y;
                    now = 1ll * now * w % mod;
                }
            }
        }
        if(sgn == -1){
            int Inv = Pow(len, mod - 2);
            for(int i = 0; i < len; i++) a[i] = 1ll * a[i] * Inv % mod;
        }
    }
}
inline int C(int x, int y){
    return 1ll * js[x] * inv[y] % mod * inv[x-y] % mod;
}
int main(){
    read(n), read(k), js[0] = inv[0] = 1;
    for(int i = 1; i <= n * k; i++){
        js[i] = 1ll * js[i-1] * i % mod;
        inv[i] = Pow(js[i], mod - 2);
    }
    f[0][0] = 1;
    for(int i = 1; i <= n; i++){
        for(int j = 0; j <= i * (k - 1); j++){
            f[i-1][j] = 1ll * f[i-1][j] * inv[j] % mod;
            if(j < k) a[j] = 1ll * inv[j] * Pow(Pow(n, mod - 2), j) % mod;
        }
        poly::timesinit(i * k);
        poly::dft(f[i-1], 1), poly::dft(a, 1);
        for(int j = 0; j < poly::len; j++)
            f[i][j] = 1ll * f[i-1][j] * a[j] % mod;
        poly::dft(f[i], -1);
        for(int j = 0; j < poly::len; j++) a[j] = 0;
        for(int j = i * (k - 1) + 1; j < poly::len; j++) f[i][j] = 0;
        for(int j = 0; j <= i * (k - 1); j++){
            f[i][j] = 1ll * f[i][j] * js[j] % mod;
            int tmp = 1ll * f[i][j] * C(n, i) % mod; 
            tmp = 1ll * tmp * Pow(n, j + 1) % mod;
            tmp = 1ll * tmp * Pow(Pow(i, mod - 2), j + 1) % mod;
            if(i & 1) up(ans, tmp); else up(ans, mod - tmp);
        }
    }
    cout << ans << endl;
    return 0;
}