LibreOJ 3315 「ZJOI2020」抽卡

继续数数。

考虑计算 i 步以后还没有连续 k 张卡的概率。
容易发现可以转化为对于所有没有连续 k 张卡的状态，计算 i 步后到达该状态的概率并求和。

设 fi 表示从 m 张卡里选择 i 张卡使得没有连续 k 张卡的方案数。
则答案即 ∞∑i=0m−1∑j=0fjj∑u=0(ju)(−1)j−u(um)i=m−1∑j=0fjj∑u=0(ju)(−1)j−u∞∑i=0(um)i=m−1∑j=0fjj∑u=0(ju)(−1)j−umm−u

后面的和式可以随便卷积一下计算，那么问题转化为计算 fi。
容易发现实际上可以将序列 a 排序后划分为若干极大连续段，对于每个极大连续段分别计算答案再分治 NTT 卷起来。

那么对于一个长度为 n 的连续段，考虑计算其中的 fi。
略加思考后不难发现这个东西就是 [xn+1](k∑j=1xj)n−i+1

即，将 n 张卡的状态视作 0/1；对于每个 0，将其与之前最长的连续 1（可以为空）合并在一起考虑；那么每一个这样的块的长度便是 1…k。
此外，需要在最后强制加入一个 0。

那么设 G(x)=x−xk+11−x，不难发现，若写出 f 的生成函数，则有 n∑i=0fizi=[xn+1]11−G(x)z

考虑对这样一个二元生成函数应用扩展拉格朗日反演，有 [xn+1]11−G(x)z=1n+1[xn]z(1−xz)2(xG−1(x))n+1

其中 G−1 为 G 的复合逆。
若求出 G−1，结合多项式求逆和多项式快速幂并展开 z(1−xz)2 便可以求得 fi。
于是问题变为求解 G−1(x)。

设 F=G−1，考虑牛顿迭代之。
设已得 G(F0(x))≡x(modxn)

在 F(x)=F0(x) 处泰勒展开 G(F(x))≡G(F0(x))+G′(F0(x))(F(x)−F0(x))(modx2n)F(x)≡F0(x)+x−G(F0(x))G′(F0(x))(modx2n)

而易得 G′(x)=kxk+1−(k+1)xk+1(1−x)2

故倍增 + 多项式全家桶即可（雾
此部分仅有一个 log，算上分治 NTT 总复杂度 O(mlog2m)。
（当然这一个 log 顶人家不知道多少次分治 NTT）

貌似有更优美做法，懒得学了（
生成函数暴力美学（

代码：

#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#define add(a,b) (a + b >= mod ? a + b - mod : a + b)
#define dec(a,b) (a < b ? a - b + mod : a - b)
using namespace std;
const int M = 2e5;
const int mod = 998244353;
inline int fpow(int a,int b)
{
    int ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
    return ret;
}
int n,m,K;
int a[M + 5];
int len[M + 5],tot;
namespace Poly
{
    const int N = 1 << 19;
    const int G = 3;
    int lg2[N + 5];
    int rev[N + 5],fac[N + 5],ifac[N + 5],inv[N + 5];
    int rt[N + 5],irt[N + 5];
    inline void init()
    {
        for(register int i = 2;i <= N;++i)
            lg2[i] = lg2[i >> 1] + 1;
        int w = fpow(G,(mod - 1) / N);
        rt[N >> 1] = 1;
        for(register int i = (N >> 1) + 1;i <= N;++i)
            rt[i] = (long long)rt[i - 1] * w % mod;
        for(register int i = (N >> 1) - 1;i;--i)
            rt[i] = rt[i << 1];
        fac[0] = 1;
        for(register int i = 1;i <= N;++i)
            fac[i] = (long long)fac[i - 1] * i % mod;
        ifac[N] = fpow(fac[N],mod - 2);
        for(register int i = N;i;--i)
            ifac[i - 1] = (long long)ifac[i] * i % mod;
        for(register int i = 1;i <= N;++i)
            inv[i] = (long long)ifac[i] * fac[i - 1] % mod;
    }
    struct poly
    {
        vector<int> a;
        inline poly(int x = 0)
        {
            x && (a.push_back(x),1);
        }
        inline poly(const vector<int> &o)
        {
            a = o,shrink();
        }
        inline poly(const poly &o)
        {
            a = o.a,shrink();
        }
        inline void shrink()
        {
            for(;!a.empty() && !a.back();a.pop_back());
        }
        inline int size() const
        {
            return a.size();
        }
        inline void resize(int x)
        {
            a.resize(x);
        }
        inline int operator[](int x) const
        {
            if(x < 0 || x >= size())
                return 0;
            return a[x];
        }
        inline void clear()
        {
            vector<int>().swap(a);
        }
        inline void ntt(int type = 1)
        {
            int n = size();
            type == -1 && (reverse(a.begin() + 1,a.end()),1);
            int lg = lg2[n] - 1;
            for(register int i = 0;i < n;++i)
                rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg),
                i < rev[i] && (swap(a[i],a[rev[i]]),1);
            for(register int w = 2,m = 1;w <= n;w <<= 1,m <<= 1)
                for(register int i = 0;i < n;i += w)
                    for(register int j = 0;j < m;++j)
                    {
                        int t = (long long)rt[m | j] * a[i | j | m] % mod;
                        a[i | j | m] = dec(a[i | j],t),a[i | j] = add(a[i | j],t);
                    }
            if(type == -1)
                for(register int i = 0;i < n;++i)
                    a[i] = (long long)a[i] * inv[n] % mod;
        }
        friend inline poly operator+(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = add(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator-(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = dec(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator*(poly a,poly b)
        {
            if(a.a.empty() || b.a.empty())
                return poly();
            int lim = 1,tot = a.size() + b.size() - 1;
            for(;lim < tot;lim <<= 1);
            a.resize(lim),b.resize(lim);
            a.ntt(),b.ntt();
            for(register int i = 0;i < lim;++i)
                a.a[i] = (long long)a[i] * b[i] % mod;
            a.ntt(-1),a.shrink();
            return a;
        }
        poly &operator+=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = add(a[i],o[i]);
            return *this;
        }
        poly &operator-=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = dec(a[i],o[i]);
            return *this;
        }
        poly &operator*=(poly o)
        {
            return (*this) = (*this) * o;
        }
        poly deriv() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() - 1);
            for(register int i = 0;i < size() - 1;++i)
                ret[i] = (long long)(i + 1) * a[i + 1] % mod;
            return poly(ret);
        }
        poly integ() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() + 1);
            for(register int i = 0;i < size();++i)
                ret[i + 1] = (long long)a[i] * inv[i + 1] % mod;
            return poly(ret);
        }
        inline poly modxn(int n) const
        {
            if(a.empty())
                return poly();
            n = min(n,size());
            return poly(vector<int>(a.begin(),a.begin() + n));
        }
        inline poly inver(int m) const
        {
            poly ret(fpow(a[0],mod - 2));
            for(register int k = 1;k < m;)
                k <<= 1,ret = (ret * (2 - modxn(k) * ret)).modxn(k);
            return ret.modxn(m);
        }
        inline poly log(int m) const
        {
            return (deriv() * inver(m)).integ().modxn(m);
        }
        inline poly exp(int m) const
        {
            poly ret(1);
            for(register int k = 1;k < m;)
                k <<= 1,ret = (ret * (1 - ret.log(k) + modxn(k))).modxn(k);
            return ret.modxn(m);
        }
        inline poly pow(int m,int k1,int k2 = -1) const
        {
            if(a.empty())
                return poly();
            if(k2 == -1)
                k2 = k1;
            int t = 0;
            for(;t < size() && !a[t];++t);
            if((long long)t * k1 >= m)
                return poly();
            poly ret;
            ret.resize(m);
            int u = fpow(a[t],mod - 2),v = fpow(a[t],k2);
            for(register int i = 0;i < m - t * k1;++i)
                ret.a[i] = (long long)operator[](i + t) * u % mod;
            ret = ret.log(m - t * k1);
            for(register int i = 0;i < ret.size();++i)
                ret.a[i] = (long long)ret[i] * k1 % mod;
            ret = ret.exp(m - t * k1),t *= k1,ret.resize(m);
            for(register int i = m - 1;i >= t;--i)
                ret.a[i] = (long long)ret[i - t] * v % mod;
            for(register int i = 0;i < t;++i)
                ret.a[i] = 0;
            return ret;
        }
    };
}
using Poly::init;
using Poly::poly;
poly f,g[M + 5],h1,h2;
int ans;
inline poly calc(int m)
{
    poly ret,powk,powk1,g,d;
    for(register int k = 1;k < m;)
        k <<= 1,
        powk1 = ((powk = ret.pow(k,K)) * ret).modxn(k),
        g.resize(2),g.a[0] = 0,g.a[1] = 1,
        g -= ((1 - ret).inver(k) * (ret - powk1)).modxn(k),
        d = (((1 - ret) * (1 - ret)).inver(k) * (powk1 * K - powk * (K + 1) + 1)).modxn(k),
        ret = (ret + g * d.inver(k)).modxn(k);
    return ret.modxn(m);
}
poly solve(int l,int r)
{
    if(l == r)
        return g[l];
    int mid = l + r >> 1;
    return solve(l,mid) * solve(mid + 1,r);
}
int main()
{
    init();
    scanf("%d%d",&m,&K),f = calc(m + 2);
    for(register int i = 0;i < f.size() - 1;++i)
        f.a[i] = f[i + 1];
    f.resize(f.size() - 1),f = f.inver(m + 1);
    for(register int i = 1;i <= m;++i)
        scanf("%d",a + i);
    sort(a + 1,a + m + 1),a[0] = -1;
    for(register int i = 1;i <= m;++i)
        a[i] == a[i - 1] + 1 ? (++len[tot]) : (len[++tot] = 1);
    for(register int i = 1;i <= tot;++i)
    {
        poly t = f.pow(len[i] + 1,len[i] + 1);
        g[i].resize(len[i] + 1);
        for(register int j = 1;j <= len[i];++j)
            g[i].a[len[i] - j + 1] = (long long)Poly::inv[len[i] + 1] * j % mod * t[len[i] - j + 1] % mod;
        g[i].a[0] = 1;
    }
    f = solve(1,tot).modxn(m),h1.resize(m),h2.resize(m);
    for(register int i = 0;i < m;++i)
        h1.a[i] = (long long)(i & 1 ? mod - 1 : 1) * Poly::ifac[i] % mod,
        h2.a[i] = (long long)m * Poly::inv[m - i] % mod * Poly::ifac[i] % mod;
    h1 = (h1 * h2).modxn(m);
    for(register int i = 0;i < m;++i)
        ans = (ans + (long long)f[i] * h1[i] % mod * Poly::fac[i] % mod) % mod;
    printf("%d\n",ans);
}