P5850 calc加强版 题解

P5850 calc加强版 题解

题意：

给定两个数 \(n,k\)，称一个长为 \(n\) 的序列 \(a\) 是合法的，当且仅当 \(a_i\in[1,k]\)，且所有 \(a_i\) 互不相同。序列 \(a\) 的权值定义为 \(\prod a_i\)。现在给定 \(m\)，对于所有 \(i\in[1,m]\)，求出 \(n=i\) 时所有合法序列的权值之和。对 \(998244353\) 取模。

\(m\le 5\times 10^5,1\le m\le k< 998244353\)

题解

首先因为序列 \(a\) 的所有位可以任意排列，所以我们最后给答案乘上 \(n!\)。现在问题就是求从 \([1,k]\) 个数中选出 \(n\) 个不同的数，求所有选法中 \(n\) 个数的乘积的和。

考虑用生成函数去刻画，就是：

\[[x^n](1+x)(1+2x)(1+3x)\ldots(1+kx) = [x^n]\prod_{i=1}^k (1+ix) \]

然后因为 \(k\) 非常大，所有求 \(k\) 个多项式的积比较困难，于是可以考虑先 \(\ln\)，然后再 \(\exp\)。

\[[x^n]\prod_{i=1}^k (1+ix) = [x^n]\exp(\sum_{i=1}^k \ln(1+ix)) \]

这里需要一个前置知识，根据麦克劳林公式，有 \(-\ln(1-x) = \sum\limits_{i=1}^\infty \frac{x^i}i\)，于是 \(\ln(1+ix) = \ln(1-(-ix))\)，代入得：

\[\begin{aligned} \sum_{i=1}^k \ln(1+ix) &= \sum_{i=1}^k\sum_{j=1}^\infty \frac{-(-ix)^j}j \\ &= \sum_{j=1}^\infty \frac{(-1)^{j+1}}j(\sum_{i=1}^ki^j)x^j \end{aligned} \]

因为我们只要前 \(n\) 项的系数，所以实际上 \(j\) 的上界只有 \(n\)，现在问题就是怎么快速计算对于每个 \(j\in[1,n]\)，\(\sum\limits_{i=1}^ki^j\) 的值，如果知道了这个，就可以做一遍多项式 \(\exp\) 解决。

考虑如何快速求解上述问题，设这个东西的 EGF 为 \(\sum\limits_{j=0}^\infty(\sum\limits_{i=0}^ki^j)\frac{x^j}{j!}\)，则：

\[\begin{aligned} \sum\limits_{j=0}^\infty(\sum\limits_{i=0}^ki^j)\frac{x^j}{j!} &= \sum_{i=0}^k\sum_{j=0}^\infty \frac{(ix)^j}{j!} \\ &= \sum_{i=0}^k e^{ix} \\ &= \frac{e^{(k+1)x}-1}{e^x-1} \end{aligned} \]

这个就可以用多项式求逆做了。于是整个题就做完了，时间复杂度 \(\mathcal{O}(n\log n)\)。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define ll long long
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
using namespace std;
const int N = (1<<20)+5,mod = 998244353;
int tr[N],n,m,lim,k;
ll f[N],g[N],inv[N],fac[N],invfac[N];
ll mi[N],tmp[N];
ll qp(ll x,int y)
{
    ll ans = 1;
    for(;y;y >>= 1,x = x*x%mod)
        if(y&1)(ans *= x) %= mod;
    return ans;
}
const int G = 3,invG = qp(3,mod-2);
int jia(int x,int y){return (x += y) >= mod?x-mod:x;}
int jian(int x,int y){return x < y?x-y+mod:x-y;}
void cpy(ll *al,ll *ar,ll *b){for(int i = 0;al != ar;al++,i++)*(b+i) = *al;}
void clr(ll *al,ll *ar,ll x = 0){for(;al != ar;al++)*al = x;}
void pri(ll *f,int n)
{
    for(int i = 0;i <= n;i++)
        printf("%lld%c",f[i],i==n?'\n':' ');
}
void init(int n)
{
    lim = 1;while(lim <= n)lim <<= 1;
    for(int i = 0;i < lim;i++)tr[i] = tr[i>>1]>>1|(i&1?lim>>1:0);
}
void NTT(ll *f,int type)
{
    for(int i = 0;i < lim;i++)if(i < tr[i])swap(f[i],f[tr[i]]);
    for(int l = 2;l <= lim;l <<= 1)
    {
        ll wn = qp(~type?G:invG,(mod-1)/l);int now = l>>1;
        for(int s = 0;s < lim;s += l)
        {
            ll w = 1;int up = s|now;
            for(int i = s;i < up;i++,w = w*wn%mod)
            {int mul = w*f[i|now]%mod;f[i|now] = jian(f[i],mul);f[i] = jia(f[i],mul);}
        }
    }
    if(~type)return ;
    int invl = qp(lim,mod-2);
    for(int i = 0;i < lim;i++)f[i] = f[i]*invl%mod;
}
ll a[N],b[N];
void times(ll *f,ll *g,ll *a,int n,int m)
{
    cpy(f,f+n+1,a);cpy(g,g+m+1,b);init(n+m);
    clr(a+n+1,a+lim);clr(b+m+1,b+lim);
    NTT(a,1);NTT(b,1);
    for(int i = 0;i <= lim;i++)a[i] = a[i]*b[i]%mod;
    NTT(a,-1);
}
void Inv(ll *f,ll *g,int n)//G = G0(2-F*G0)
{
    if(!n)return (void)(g[0] = qp(f[0],mod-2));
    Inv(f,g,n>>1);init(n << 1);
    cpy(f,f+n+1,a);clr(a+n+1,a+lim);clr(g+n+1,g+lim);
    NTT(a,1);NTT(g,1);
    for(int i = 0;i < lim;i++)
        g[i] = g[i]*(2-a[i]*g[i]%mod+mod)%mod;
    NTT(g,-1);clr(g+n+1,g+lim);
}
void dao(ll *f,ll *g,int n)
{for(int i = 1;i <= n;i++)g[i-1] = f[i]*i%mod;g[n] = 0;}
void jifen(ll *f,ll *g,int n)
{for(int i = n;~i;i--)g[i+1] = f[i]*inv[i+1]%mod;g[0] = 0;}
ll invf[N];
void ln(ll *f,ll *g,int n)//G = jifen(F'[i]/F[i])
{
    Inv(f,invf,n);dao(f,g,n); 
    times(g,invf,g,n,n);clr(g+n+1,g+lim);
    jifen(g,g,n-1);
}
ll lng[N];
void exp(ll *f,ll *g,int n)//G = G0(1-lnG0+F)
{
    if(!n)return (void)(g[0] = 1);
    exp(f,g,n>>1);ln(g,lng,n);
    for(int i = 0;i <= n;i++)lng[i] = jian(!i+f[i],lng[i]);
    times(g,lng,g,n,n);clr(g+n+1,g+lim);
}
char buf[1<<21],*p1,*p2;
inline int rd()
{
    char c;int f = 1;
    while(!isdigit(c = getchar()))if(c=='-')f = -1;
    int x = c-'0';
    while(isdigit(c = getchar()))x = x*10+(c^48);
    return x*f;
}
int main()
{
    fac[0] = invfac[0] = inv[0] = inv[1] = 1;
    for(int i = 2;i < N;i++)inv[i] = inv[mod%i]*(mod-mod/i)%mod;
    for(int i = 1;i < N;i++)fac[i] = fac[i-1]*i%mod,invfac[i] = invfac[i-1]*inv[i]%mod;
    k = rd();n = rd();
    for(int i = 0;i <= n;i++)mi[i] = invfac[i+1];
    Inv(mi,tmp,n);
    for(int i = 0;i <= n;i++)mi[i] = qp(k+1,i+1)*invfac[i+1]%mod;
    times(mi,tmp,mi,n,n);
    for(int i = 1;i <= n;i++)f[i] = inv[i]*(i&1?1:mod-1)%mod*mi[i]%mod*fac[i]%mod;
    exp(f,g,n);
    for(int i = 1;i <= n;i++)printf("%lld\n",g[i]*fac[i]%mod);
    return 0;
}