Luogu P7445「EZEC-7」线段树
显然一个点是否被( ext{push_down})仅取决于所有完全包含它的操作区间权值之和
那么可以考虑对于每个节点计算概率，然后累加
反向计算一个节点不被( ext{push_down})的概率，即权值之和为(0)的概率
而每个节点有自己被覆盖的概率，即(p_i=cfrac{lcdot (n-r+1)}{n(n+1)/2})
而覆盖的次数(c)决定了这个概率贡献的权值，即(p_i^c(1-p_i)^{m-c})
由此得到一个思路：
先通过计算得到(k)次覆盖权值为0的函数(A(x))
容易发现这样得到每个点的概率为：(A(cfrac{p_i}{1-p_i})(1-p_i)^m)
Naive地带入多点求值，就能暴力得到
计算(k)次被覆盖权值和为0的方案数
容易发现就是
(displaystyle [x^0](sum_{i=-1}^V x^i)^k=[x^0](x^{-1}frac{1-x^{V+2}}{1-x})^k)
暴力展开这个式子会需要对于((x^{V+2}-1)^k)有用的项只有(frac{k}{V+2})项
即原式(displaystyle =[x^k](frac{1-x^{V+2}}{1-x})^k=sum_{i=0}^{k}inom{2k-i-1}{k-1}[x^i](1-x^{V+2})^k)
(displaystyle =sum_{i=0}^{frac{k}{V+2}} inom{2k-i(V+2)-1)}{k-1} (-1)^i inom{k}{i})
(第一个组合数是组合意义插板，第二个是二项展开)
求(k)的权值需要(O(frac{k}{V}))，并不好直接卷积优化
由于涉及了类似([x^n]G^k(x))的形式，考虑用 "另类拉格朗日反演" 求解
处理一下(x)的负指数，设(displaystyle F(x)=sum _{i=0}^{V+1}x^i)，转化为([x^0](frac{F(x)}{x})^k)
然而不管是(F(x))还是(frac{F(x)}{x})都不存在复合逆，但是(frac{x}{F(x)})有
设(G(x))为(frac{x}{F(x)})的复合逆
(displaystyle [x^0](frac{F(x)}{x})^k=[x^0](frac{x}{F(x)})^{-k}=[x^{k}]frac{xG'(x)}{G(x)})
求解(G(x))即可得到所有(k)的值
(G(x))为(cfrac{x}{F(x)}=cfrac{x (x-1)}{x^{V+2}-1})的复合逆
即满足(cfrac{G(G-1)}{G^{V+2}-1}=x)
(xG^{V+2}-G^2+G-x=0)
这个形式还是比较易于进行牛顿迭代的
(f(z)=xz^{V+2}-z^2+z-x)
(f'(z)=(V+2)xz^{V+1}-2z+1)
(displaystyle A=B-frac{f(B)}{f'(B)}=B-frac{xB^{V+2}-B^2+B-x}{(V+2)xB^{V+1}-2B+1})
边界条件为([x^0]G(x)=0)，([x^1]G(x)=1)
结果牛顿迭代需要多项式快速幂
有关多点求值的优化
由于我们并不需要知道每个点被操作的概率，只需要一个求和，因此可以对于每一项求出
设(a_i=cfrac{p_i}{1-p_i},b_i=(1-p_i)^m)，容易发现实际上求出的式子是
(A(cfrac{p_i}{1-p_i})(1-p_i)^m=sum A_jsum_i a_i^j b_i)
对于每个(j)求解，就是
(displaystyle [x^j]sum _i frac{b_i}{1-a_ix})
可以分治( ext{NTT})通分，也就是写成下式
(displaystyle frac{1}{prod (1-a_ix)} sum b_i prod_{i!=j} (1-a_jx))
右边就是一个经典的分治( ext{NTT})问题，再加上一次求逆即可
好像也不一定快吧
接下来就是套板板时间
const int N=1<<19,P=998244353;
int n,m,k;
ll qpow(ll x,ll k=P-2) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
typedef vector <int> V;
int rev[N],w[N];
int Inv[N+1],I[N+1],J[N+1];
void Init_w() {
int N=1; while(N<=max(n,m+1)*2+4) N<<=1;
int t=qpow(3,(P-1)/N);
w[N/2]=1;
rep(i,N/2+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,N/2-1,1) w[i]=w[i<<1];
Inv[0]=Inv[1]=1;
rep(i,2,N) Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
rep(i,*I=*J=1,N) {
I[i]=1ll*I[i-1]*Inv[i]%P;
J[i]=1ll*J[i-1]*i%P;
}
}
int Init(int n){
int R=1,c=-1;
while(R<=n) R<<=1,c++;
rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
return R;
}
/*
被隐藏的部分:!!
NTT
operator *
operator +
operator -
*/
V operator ~ (V a) {
int n=a.size();
if(n==1) return V{(int)qpow(a[0],P-2)};
V b=a; b.resize((n+1)/2); b=~b;
int R=Init(n*2);
NTT(R,a,1),NTT(R,b,1);
rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
NTT(R,a,-1),a.resize(n);
return a;
}
void Exp_Solve(V &A,V &B,int l,int r){
static int X[N],Y[N];
if(l==r) {
B[l]=1ll*B[l]*Inv[l]%P;
return;
}
int mid=(l+r)>>1;
Exp_Solve(A,B,l,mid);
int R=Init(r-l+2);
rep(i,0,R) X[i]=Y[i]=0;
rep(i,l,mid) X[i-l]=B[i];
rep(i,0,r-l-1) Y[i+1]=A[i];
NTT(R,X,1),NTT(R,Y,1);
rep(i,0,R-1) X[i]=1ll*X[i]*Y[i]%P;
NTT(R,X,-1);
rep(i,mid+1,r) B[i]+=X[i-l],Mod1(B[i]);
Exp_Solve(A,B,mid+1,r);
}
V Deri(V a){
rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
a.pop_back();
return a;
}
V Integ(V a) {
a.pb(0);
drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Inv[i]%P;
return a[0]=0,a;
}
V operator << (V A,const int &x){
A.resize(A.size()+x);
drep(i,A.size()-1,x) A[i]=A[i-x];
rep(i,0,x-1) A[i]=0;
return A;
}
V operator >> (V A,const int &x){
rep(i,x,A.size()-1) A[i-x]=A[i];
A.resize(A.size()-x);
return A;
}
V Ln(V a){
int n=a.size();
a=Deri(a)*~a,a.resize(n-1);
return Integ(a);
}
V Exp(V F){
int n=F.size(); F=Deri(F);
V A(n); A[0]=1;
Exp_Solve(F,A,0,n-1);
return A;
}
V Pow(V x,int k) {
int d=0,n=x.size();
while(d<n && !x[d]) d++;
if(1ll*d*k>=n){
rep(i,0,x.size()-1) x[i]=0;
return x;
}
x=x>>d,x.resize(n),x=Ln(x);
rep(i,0,n-1) x[i]=1ll*x[i]*k%P;
x=Exp(x)<<(d*k),x.resize(n);
return x;
}
V Evaluate(V F,V X){
static int ls[N<<1],rs[N<<1],cnt;
static V T[N<<1];
static auto TMul=[&] (V F,V G){
int n=F.size(),m=G.size();
reverse(G.begin(),G.end());
int R=Init(n);
NTT(R,F,1),NTT(R,G,1);
rep(i,0,R-1) F[i]=1ll*F[i]*G[i]%P;
NTT(R,F,-1); V T(n-m+1);
rep(i,0,n-m) T[i]=F[i+m-1];
return T;
};
static function <int(int,int)> Build=[&](int l,int r) {
int u=++cnt; ls[u]=rs[u]=0;
if(l==r) {
T[u]=V{1,P-X[l]};
return u;
}
int mid=(l+r)>>1;
ls[u]=Build(l,mid),rs[u]=Build(mid+1,r);
T[u]=T[ls[u]]*T[rs[u]];
return u;
};
int n=F.size(),m=X.size();
cmax(n,m),F.resize(n),X.resize(n);
cnt=0,Build(0,n-1);
F.resize(n*2+1),T[1]=TMul(F,~T[1]);
int p=0;
rep(i,1,cnt) if(ls[i]) {
swap(T[ls[i]],T[rs[i]]);
int R=Init(T[i].size()),n=T[i].size(),m1=T[ls[i]].size(),m2=T[rs[i]].size();
NTT(R,T[i],1);
reverse(T[ls[i]].begin(),T[ls[i]].end()); reverse(T[rs[i]].begin(),T[rs[i]].end());
NTT(R,T[ls[i]],1); NTT(R,T[rs[i]],1);
rep(j,0,R-1) {
T[ls[i]][j]=1ll*T[ls[i]][j]*T[i][j]%P;
T[rs[i]][j]=1ll*T[rs[i]][j]*T[i][j]%P;
}
NTT(R,T[ls[i]],-1); NTT(R,T[rs[i]],-1);
rep(j,0,n-m1) T[ls[i]][j]=T[ls[i]][j+m1-1];
T[ls[i]].resize(n-m1+1);
rep(j,0,n-m2) T[rs[i]][j]=T[rs[i]][j+m2-1];
T[rs[i]].resize(n-m2+1);
} else X[p++]=T[i][0];
X.resize(m);
return X;
}
V operator * (V A,const int &x){
rep(i,0,A.size()-1) A[i]=1ll*A[i]*x%P;
return A;
}
V Newton(int n){
if(n==1) return V{0,1};
V G=Newton((n+1)/2); G.resize(n);
V T=Pow(G,k+1);
V A=((G*T)<<1)-G*G+G-V{0,1},B=(T<<1)*(k+2)-G*2+V{1};
A.resize(n+1),B.resize(n+1),A=A*~B,A.resize(n+1);
return G-A;
}
V X,Y;
void Build(int l,int r){
if(l==r) return;
int prob=1ll*l*(n-r+1)%P*Inv[n]%P*Inv[n+1]%P*2%P;
Y.pb(P+1-prob);
X.pb(prob*qpow(P+1-prob)%P);
int mid=(l+r)>>1;
Build(l,mid),Build(mid+1,r);
}
int main(){
n=rd(),m=rd(),k=rd(),Init_w();
V F=Newton(m+1);
F=Deri(F)*~(F>>1),F.resize(m+1);
int t=1,inv=qpow(k+2);
rep(i,0,m) {
F[i]=1ll*F[i]*t%P*J[m]%P*I[i]%P*I[m-i]%P;
t=1ll*t*inv%P;
}
Build(1,n);
X=Evaluate(F,X);
int ans=0;
rep(i,0,X.size()-1) ans=(ans+P+1-X[i]*qpow(Y[i],m))%P;
Mod2(ans),printf("%d
",ans);
}