P6667 [清华集训2016] 如何优雅地求和 -Binomial Sum

发布时间 2023-09-23 11:15:50作者: British_Union

题面

有一个多项式函数 \(f(x)\),最高次幂为 \(x^m\),定义变换 \(Q\)

\[Q(f,n,x)=\sum_{k=0}^{n}f(k)\binom{n}{k}x^k(1-x)^{n-k} \]

现在给定函数 \(f\)\(n,x\),求 \(Q(f,n,x)\bmod998244353\)
出于某种原因,函数 \(f\) 由点值形式给出,即给定 \(a_0,a_1,⋯,a_m\)\(m+1\) 个数,\(f(x)=a_x\)。可以证明该函数唯一。

Solution

不知道为什么没有 Binomial Sum 的题解?
Binomial Sum 算法要求已求得 \(\displaystyle\sum_{i=0}^m a_i[x^i]G(x)^k,k\in[1,n]\)
然后可以求得 \(\displaystyle\sum_{i=0}^m a_i[x^i]F(G(x))\)
构造 \(G(x)=e^x,F(x)=(ax+1-a)^n\)。(避免混淆,这里把 \(x\) 写作 \(a\),但是注意 \(a_i\)\(a\) 的区别)
接下来说明这个东西的正确性:
\([x^j]F(G(x))=(ae^x+1-a)^n =\displaystyle\sum_{i=0}^n \binom{n}{i}a^i(1-a)^{n-i}[x^j]e^{ix}\)
\([x^j]e^{ix}=\dfrac{i^j}{j!}\),所以原式\(=\displaystyle\sum_{i=0}^n \binom{n}{i}a^i(1-a)^{n-i}\dfrac{i^j}{j!}\)
不妨设 \([x^j]f(x)=\frac{a_j}{j!}\),那么 \(\displaystyle \sum_{j=0}^ma_j[x^j]G(x)^k=f(k)\)。(注意 \(a_i\) 并非给定)
式子 \(\displaystyle\sum_{j=0}^ma_j[x^j]\sum_{i=0}^n \binom{n}{i}a^i(1-a)^{n-i}\dfrac{i^j}{j!}\)
\(=Q(f,n,x)\)
然后按照 Binomial Sum 的现成方法来就可以了。
这里容易发现 \(F\) 满足的一个微分方程:
\((ax+1-a)F'(x)-anF(x)=0\)

\(G(0)=1\),所以设 \(\mathcal{F}(x)=F(x)\bmod(x-1)^{m+1}\),或\(\mathcal{F}(x+1)=F(x+1)\bmod x^{m+1}=(ax+1)^n \bmod x^{m+1}\)

微分方程变换为 \((ax+1)F'(x+1)-anF(x+1)=0\),以预备换为关于 \(\mathcal{F}\) 的方程。

考察 \((ax+1)\mathcal{F}'(x+1)-an\mathcal{F}(x+1)\),其不合法的项从哪里来?

显然只有 \(\mathcal{F}'\) 从其原函数的 \(m+1\) 次跨向 \(m\) 次中来。

所以 \((ax+1)\mathcal{F}'(x+1)\),具体地说是后面那个 \(1\mathcal{F}'(x+1)\) 获得了贡献。但是实际上不存在。所以要在右边减去。

所以右边等于 \(-x^m[x^m]F'(x+1)=-na^{m+1}\binom{n-1}{m}x^m\)

于是右边是一个可以变化、展开的项,令 \(x-1\to x\)

然后得到 \((ax+1-a)\mathcal{F}'(x)-an\mathcal{F}(x)=-na^{m+1}\binom{n-1}{m}(x-1)^{m}\)

于是提取系数,具体来说:

\([x^i]\mathcal{F}(x)=f_i\)\([x^i]\) 左式 \(=(1-a)(i+1)f_{i+1}+aif_i-anf_i\)

同时 \([x^i]\) 右式 \(=-na^{m+1}\binom{n-1}{m}\binom{m}{i}(-1)^{m-i}\)

由于两边多项式相等,系数也相等,左式等于右式,得出 \(f_i\) 递推式。

此时计算得到 \(f_0\) 的值是不明智的。注意到 \(f_{m+1}=0\),倒着推即可。注意特判 \(n\le m\)

精细实现(预处理连续逆元)可以做到 \(O(m)\),不精细实现则是小常数 \(O(m\log m)\),吊打 NTT。

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn=2e4+5,mod=998244353;
int sqr(int x){
	return x*x%mod;
}
int qp(int a,int b){
	if(b==0)return 1;
	return b==1?a:(sqr(qp(a,b>>1))*(b%2?a:1)%mod);
}
int n,m,a,g[maxn],fac[maxn],ifac[maxn];
int f[maxn];
int pow1(int x){
	return x%2?mod-1:1;
}
int Mod(int x){
	return (x%mod+mod)%mod;
}
int C(int a,int b){
	return fac[a]*ifac[b]%mod*ifac[a-b]%mod;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	cin>>n>>m>>a;
	for(int i=0;i<=m;i++)cin>>g[i];
	fac[0]=1;
	for(int i=1;i<=m;i++)fac[i]=fac[i-1]*i%mod;
	ifac[m]=qp(fac[m],mod-2);
	for(int i=m-1;i>=0;i--)ifac[i]=ifac[i+1]*(i+1)%mod;
	int res=0;
	if(n<=m){
		for(int i=0;i<=n;i++){
			res+=g[i]*C(n,i)%mod*qp(a,i)%mod*qp(mod+1-a,n-i)%mod;
			res%=mod;
		}
		cout<<res<<endl;
		return 0;
	}
	int bn=ifac[m],am=qp(a,m+1);//(n-1,m),a^m
	for(int i=1;i<=m;i++){
		bn=bn*(n-i)%mod;
	}
	f[m+1]=0;
	for(int i=m;i>=0;i--){
		f[i]=qp(a*(n-i)%mod,mod-2)*Mod((i+1)*(mod+1-a)%mod*f[i+1]%mod-n*am%mod*bn%mod*C(m,i)%mod*pow1(m-i+1)%mod)%mod;
	}
	for(int i=0;i<=m;i++){
		res+=g[i]*f[i]%mod;
		res%=mod;
	}
	cout<<Mod(res)<<endl;
	return 0;
}