洛谷 P4463 - [集训队互测 2012] calc（多项式）
竟然能独立做出 jxd 互测的题（及其加强版），震撼震撼（（（故写题解以祭之
首先由于 \(a_1,a_2,\cdots,a_n\) 互不相同，故可以考虑求出所有集合 \(S=\{a_1,a_2,\cdots,a_n\}\) 的权值之和，然后答案乘上 \(n!\)。
那么怎么求这个权值之和呢？首先有一个非常 naive 的 DP，\(dp_{i,j}\) 表示 \(1\sim i\) 中选了 \(j\) 个数，可得的集合的权值之和，那么显然有 \(dp_{i,j}=dp_{i-1,j}+i\times dp_{i-1,j-1}\)，最终答案即为 \(dp_{k,n}\)。
不难发现这个 DP 可以写成多项式的形式，记 \(F_i(x)\) 为 \(dp_i\) 的 OGF，那么有 \(F_i(x)=F_{i-1}(x)\times (1+ix)\)，即 \(F_i(x)=\prod\limits_{j=1}^i(1+jx)\)，也就是说最终答案即为 \(n![x^n]\prod\limits_{j=1}^k(1+jx)\)。
那么怎么计算这个东西呢？我的第一想法是分治计算上式的值，但很快否决掉了这个想法，因为此题的 \(k\) 高达 \(10^9\)，根本无法分治计算。注意到像 \(1+jx\) 这类非常简洁的式子，它的 \(\ln\) 也是可以以非常简洁的形式表示出来的，为 \(jx-\dfrac{(jx)^2}{2}+\dfrac{(jx)^3}{3}-\cdots\)，故考虑借鉴 P4389 付公主的背包 的套路，求个 \(\ln\) 加起来再 \(\exp\) 回去，即 \(\exp(\sum\limits_{j=1}^k\ln(1+jx))=\exp(\sum\limits_{j=1}^k\sum\limits_{i}(-1)^{i+1}\dfrac{(jx)^i}{i})\)，将 \(\exp\) 里面的东西简单化个简即可得到 \(\sum\limits_{i}\dfrac{(-1)^{i+1}\sum\limits_{j=1}^nj^i}{i}x^i\)，发现瓶颈在于自然数 \(i\) 次方和，这个显然可以通过 lagrange 插值求出。复杂度 \(n^2\)。
最后此题（原题）的模数比较恶心，但是由于数据范围比较小，暴力卷积就能过，于是就懒得写 MTT 了（实际上是懒癌症晚期（（（
btw 由于以数组形式实现的多项式非常容易出现全局变量不清空的问题，因此从这道题开始我将改用 vector 实现多项式全家桶，实测非常好写（bushi）并且效率也不低。
const int MAXN=500;
const int MAXP=1024;
int k,n,mod,inv[MAXP+5];
int qpow(int x,int e=mod-2){
int ret=1;
for(;e;e>>=1,x=1ll*x*x%mod) if(e&1) ret=1ll*ret*x%mod;
return ret;
}
vector<int> conv(vector<int> a,vector<int> b){
vector<int> c(a.size()+b.size(),0);
for(int i=0;i<a.size();i++) for(int j=0;j<b.size();j++)
c[i+j]=(c[i+j]+1ll*a[i]*b[j]%mod)%mod;
return c;
}
vector<int> getinv(vector<int> a,int n){
vector<int> b(n,0);b[0]=qpow(a[0]);
for(int i=2;i<=n;i<<=1){
vector<int> c(a.begin(),a.begin()+i);
vector<int> d(b.begin(),b.begin()+(i>>1));
d=conv(d,d);c=conv(c,d);
for(int j=0;j<i;j++) b[j]=(2*b[j]%mod-c[j]+mod)%mod;
} return b;
}
vector<int> direv(vector<int> a,int n){
vector<int> b(n,0);
for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;
return b;
}
vector<int> inter(vector<int> a,int n){
vector<int> b(n,0);
for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*inv[i]%mod;
return b;
}
vector<int> getln(vector<int> a,int n){
vector<int> a_=direv(a,n),b=getinv(a,n);
b=conv(a_,b);b=inter(b,n);return b;
}
vector<int> getexp(vector<int> a,int n){
vector<int> b(n,0);b[0]=1;
for(int i=2;i<=n;i<<=1){
vector<int> c(b.begin(),b.begin()+i);
vector<int> d(b.begin(),b.begin()+(i>>1));
c=getln(c,i);for(int j=0;j<i;j++) c[j]=(a[j]-c[j]+mod)%mod;
(c[0]+=1)%=mod;d=conv(d,c);
for(int j=0;j<i;j++) b[j]=d[j];
} return b;
}
int ff[MAXN+5],pre[MAXN+5],suf[MAXN+5],fac[MAXN+5],ifac[MAXN+5];
void init_fac(int n){
fac[0]=ifac[0]=ifac[1]=1;
for(int i=2;i<=n;i++) ifac[i]=1ll*ifac[mod%i]*(mod-mod/i)%mod;
for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%mod,ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;
}
int getsum(int t,int k){
int ans=0;
ff[0]=0;for(int i=1;i<=k+1;i++) ff[i]=(ff[i-1]+qpow(i,k))%mod;
pre[0]=t;for(int i=1;i<=k+1;i++) pre[i]=1ll*pre[i-1]*(t-i+mod)%mod;
suf[k+2]=1;for(int i=k+1;~i;i--) suf[i]=1ll*suf[i+1]*(t-i+mod)%mod;
for(int i=0;i<=k+1;i++){
int mul=1;
if(i!=0) mul=1ll*mul*pre[i-1]%mod;
if(i!=k+1) mul=1ll*mul*suf[i+1]%mod;
mul=1ll*mul*ifac[i]%mod;mul=1ll*mul*ifac[k+1-i]%mod;
if((k+1-i)&1) mul=(mod-mul)%mod;
// printf("%d %d\n",i,mul);
ans=(ans+1ll*mul*ff[i])%mod;
} return ans;
}
int main(){
scanf("%d%d%d",&k,&n,&mod);init_fac(n+1);
// printf("%d\n",getsum(10,4));
int LEN=1;while(LEN<=n) LEN<<=1;
for(int i=1;i<=LEN;i++) inv[i]=qpow(i);
vector<int> a(LEN);
for(int i=1;i<=n;i++){
if(i&1) a[i]=1ll*getsum(k,i)*inv[i]%mod;
else a[i]=(mod-1ll*getsum(k,i)*inv[i]%mod)%mod;
// printf("%d %d\n",a[i],getsum(k,i));
} a=getexp(a,LEN);int ret=a[n];
for(int i=1;i<=n;i++) ret=1ll*ret*i%mod;
printf("%d\n",ret);
return 0;
}
/*
3 2 10007
*/
/*
9 7 10007
*/
在原题中 \(n\) 的数据范围只有 \(500\)，故暴力插值也能过，但如果 \(n\) 扩大到了 \(5\times 10^5\) 怎么办呢？
发现在上述算法中我们使用的是拉格朗日插值法求自然数 \(i\) 次幂的和，并且对所有 \(i\in[1,n]\) 都插了一遍，效率低下。那么有什么办法能够快速对 \(i\in[1,n]\) 求出 \(\sum\limits_{j=1}^kj^i\) 呢？
设序列 \(b_i=\sum\limits_{j=1}^kj^i\)，考虑 \(b\) 序列的 EGF，它等于 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^k\dfrac{j^ix^i}{i!}=\sum\limits_{j=1}^k\sum\limits_{i=1}^n\dfrac{(jx)^i}{i!}=\sum\limits_{j=1}^ke^{jx}=\dfrac{e^{(k+1)x}-1}{e^x-1}\)，噫，好，可以直接算了。多项式求逆求出上述式子的值即可。虽然底下分母 \(e^x-1\) 的常数项为 \(0\)，不存在乘法逆，但是 \(e^{(k+1)x}-1\) 的常数项也为 \(0\)，故可以分子分母同除 \(x\)，这样就可用常规的多项式求逆计算了。
注意点：多项式 exp 时必须保证常数项为 \(0\)，否则会出现奇怪的错误。
const int pr=3;
const int MOD=998244353;
const int ipr=(MOD+1)/3;
const int MAXP=1<<20;
int qpow(int x,int e=MOD-2){
int ret=1;
for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
return ret;
}
int n,m,rev[MAXP+5],inv[MAXP+5],fac[MAXP+5],ifac[MAXP+5];
void init_fac(int n){
fac[0]=ifac[0]=ifac[1]=1;
for(int i=2;i<=n;i++) ifac[i]=1ll*ifac[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%MOD,ifac[i]=1ll*ifac[i-1]*ifac[i]%MOD;
}
void NTT(vector<int> &a,int len,int type){
int lg=31-__builtin_clz(len);
for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=2;i<=len;i<<=1){
int W=qpow((type<0)?ipr:pr,(MOD-1)/i);
for(int j=0;j<len;j+=i){
int w=1;
for(int k=0;k<(i>>1);k++,w=1ll*w*W%MOD){
int X=a[j+k],Y=1ll*a[(i>>1)+j+k]*w%MOD;
a[j+k]=(X+Y)%MOD;a[(i>>1)+j+k]=(X-Y+MOD)%MOD;
}
}
}
if(type<0) for(int i=0;i<len;i++) a[i]=1ll*a[i]*inv[len]%MOD;
}
vector<int> conv(vector<int> a,vector<int> b){
int len=1;while(len<a.size()+b.size()) len<<=1;
a.resize(len,0);b.resize(len,0);NTT(a,len,1);NTT(b,len,1);
for(int i=0;i<a.size();i++) a[i]=1ll*a[i]*b[i]%MOD;
NTT(a,len,-1);return a;
}
vector<int> getinv(vector<int> a,int n){
vector<int> b(n,0);b[0]=qpow(a[0]);
for(int i=2;i<=n;i<<=1){
vector<int> c(a.begin(),a.begin()+i);
vector<int> d(b.begin(),b.begin()+(i>>1));
d=conv(d,d);c=conv(c,d);
for(int j=0;j<i;j++) b[j]=(2*b[j]%MOD-c[j]+MOD)%MOD;
} return b;
}
vector<int> direv(vector<int> a,int n){
vector<int> b(n,0);
for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%MOD;
return b;
}
vector<int> inter(vector<int> a,int n){
vector<int> b(n,0);
for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*inv[i]%MOD;
return b;
}
vector<int> getln(vector<int> a,int n){
vector<int> a_=direv(a,n),b=getinv(a,n);
b=conv(a_,b);b=inter(b,n);return b;
}
vector<int> getexp(vector<int> a,int n){
vector<int> b(n,0);b[0]=1;
for(int i=2;i<=n;i<<=1){
vector<int> c(b.begin(),b.begin()+i);
vector<int> d(b.begin(),b.begin()+(i>>1));
c=getln(c,i);for(int j=0;j<i;j++) c[j]=(a[j]-c[j]+MOD)%MOD;
(c[0]+=1)%=MOD;d=conv(d,c);
for(int j=0;j<i;j++) b[j]=d[j];
} return b;
}
int main(){
scanf("%d%d",&n,&m);int LEN=1;
while(LEN<=m+1) LEN<<=1;init_fac(LEN);
for(int i=1;i<=LEN<<1;i++) inv[i]=qpow(i);
vector<int> a(LEN,0),b(LEN,0);
for(int i=1,mul=n+1;i<=m+1;i++,mul=1ll*mul*(n+1)%MOD){
a[i-1]=ifac[i];b[i-1]=1ll*mul*ifac[i]%MOD;
} a=getinv(a,LEN);a=conv(a,b);
for(int i=1;i<=m+1;i++) a[i]=1ll*a[i]*fac[i-1]%MOD;
for(int i=2;i<=m+1;i+=2) a[i]=(MOD-a[i])%MOD;
a[0]=0;a=getexp(a,LEN);
for(int i=1;i<=m;i++) printf("%d\n",1ll*a[i]*fac[i]%MOD);
return 0;
}
