LOJ3075 「2019 集训队互测 Day 3」组合数求和
题意：
令 \(f(j)=\sum\limits_{i=0}^{n-1}\dbinom{id}{j}\) ，对于 \(0\le j <m\) ，分别求出 \(f(j)\) 。答案对 \(M\) 取模。
\(10^8\le M\le 10^9\) (不保证是素数) 。
\(1\le nd\le 10^9,1\le md\le 3*10^6,1\le d\le 100\) 。
做法：
step1
首先答案能用生成函数的角度求出：令 \(C(x)=\sum\limits_{i=0}^{n-1} (1+x)^{id}\)，则 \(f(j)=[x^j]C(x)\) 。 就是求 \(C(x)\) 的前 \(m\) 项系数。
考虑这是等比数列求和，令 \(A(x)=\frac{(1+x)^{nd}-1}{x},B(x)=\frac{(1+x)^{d}-1}{x}\) ，则 \(C(x)=\frac{A(x)}{B(x)}\) 。
step2
\(A\) 和 \(B\) 都是容易求出的，把模数拆成 \(M=\prod p_i^{c_i}\) 的形式，求出模 \(p_i^{c_i}\) 时的一行组合数是不难的（记录一下 \(p_i\) 的次数，其他的正常乘），然后再 CRT 合并。
这一部分复杂度大概是 \(O((m+d)\log n)\) 的。
step3
令 \(f_i=[x^i]A(x),g_i=[x^i]B(x),h_i=[x^i]C(x)\) 。
先考虑 \(gcd(d,M)=1\) 这档分， \(g_0h_i=f_i-\sum\limits_{j=1}^{d-1}g_jh_{i-j}\) 。
因为 \(g_0=d\) ，所以可以求逆元，我们就能 \(O(md)\) 算出来。
\(gcd(d,M)\ne 1\) 怎么办？到这里我就不会做了/kk
step4
还是把 \(M\) 拆成 \(\prod p_i^{c_i}\) 的形式，把模 \(p^C\) 的答案算出来，再 CRT 合并。注意 \(d\) 不是 \(p_i\) 倍数的部分，可以乘在一起，用上面的方法处理。
对于 \(d\) 是 \(p\) 的倍数，找出最小的 \(t\) 使 \(g_t\) 不是 \(p\) 的倍数。由于 \(g_{d-1}=1\) ，我们是一定能找出来的。
则 \(h_ig_t=f_{i+t}-\sum\limits_{j\ne t}g_jh_{i+t-j}\) 。(称为式子 \(P\))
\(g_t\) 是能求逆元的，但是我们的式子中存在 “高” 贡献给 "低" 的情况，这是不好的。
我们以 \(i\) 为主元，维护 \(a\) 和 \(c\)，表示 \(h_i=\sum a_j f_{i+j}+\sum c_j h_{i+j}\) 。
我们把式子 \(P\) 用上面的方式表达出来，记为式子 \(Q\) 。
接下来展示一种方法，使得能把一个初始 \(=Q\) 的式子 \(R\) ，通过变换，使得 \(R\) 中对于 \(j>0\) 都有 \(c_j\) 为 \(p^C\) 倍数。利用 \(R\) ，就能以 \(O(m|R|)\) 的时间递推出 \(h\) 。
(只存在 “低” 贡献给 “高” 了)
根据 \(t\) 的最小化的定义，变换前 \(R\) 满足：对于 \(j>0\) 都有 \(c_j\) 为 \(p^1\) 的倍数。
如果对于第 \(i\) 轮操作，初始对于 \(j>0\) 都有 \(c_j\) 为 \(p^i\) 倍数 ,我们能通过变换使得 \(c_j\) 都为 \(p^{i+1}\) 的倍数。则通过 \(C-1\) 轮操作使 \(R\) 合法。
而操作的方法就是，从大到小枚举 \(j>0\) ，如果当前 \(c_j\) 不是 \(p^{i+1}\) 的倍数，就把 \(h_{i+j}\) 利用 \(Q\) 展开。（展开 \(h_{i+j}\) ，也就是把 \(Q\) 中的 \(i\) 用 \(i+j\) 替换）
此时，\(k\geq j\) 都满足是 \(p^{i+1}\) 的倍数，且 \(k<j\) 仍满足是 \(p^i\) 的倍数。
则每一轮操作的时间复杂度是 \(O(d|R|)\) 。因为每次 \(|R|\) 会增加至多 \(d-1\) ，所以求 \(R\) 总复杂度是 \(O(C^2d^2)\) ，最后 \(|R|\) 能达到 \(O(Cd)\) 的长度。
得出模 \(p^C\) 的复杂度是 \(O(Cmd+C^2d^2)\) 。
最后得到的 \(R\) 有可能 \(c_0\ne 0\) ，但它一定是 \(p\) 的倍数，所以 \(gcd(1-c_0,p^C)=1\) ，再求一下逆元就好了。
总结
总复杂度为 \(O(md\log M+d^2\log^2 M+(m+d)\log n)\) 。
自己没想到的东西：
-
拆模数，求模数为 \(p_i^{c_i}\) 的答案后 CRT 合并。这是很常见的方法。
-
处理不好求逆元的多项式除法的神秘方法。
代码：
#include<bits/stdc++.h>
#define ll long long
using namespace std;
int n,m,d,M,mod;
int P[10],C[10],t;
void exgcd(int a,int b,ll &x,ll &y){
if(!b){x=1,y=0;return;}
exgcd(b,a%b,y,x),y-=(a/b)*x;
}
int inv(int a,int b){ll x,y;exgcd(a,b,x,y);return (x%b+b)%b;}
int A[3010000],B[110];
void so(int p,int c,int xs){
int mo=1;
for(int i=1;i<=c;i++)mo*=p;
int Z=1,PX=0;
for(int i=1;i<=m+10010;i++){
int a=n-i+1;
if(!a)break;
while(!(a%p))a/=p,PX++;
Z=1ll*Z*a%mo;
a=i;while(!(a%p))a/=p,PX--;
Z=1ll*Z*inv(a,mo)%mo;
int az=Z,pt=PX;
if(pt>=c)az=0;
else{while(pt--)az=1ll*p*az%mod;}
(A[i-1]+=1ll*az*xs%mod)%=mod;
}
Z=1,PX=0;
for(int i=1;i<=d;i++){
int a=d-i+1;
while(!(a%p))a/=p,PX++;
Z=1ll*Z*a%mo;
a=i;while(!(a%p))a/=p,PX--;
Z=1ll*Z*inv(a,mo)%mo;
int az=Z,pt=PX;
if(pt>=c)az=0;
else{while(pt--)az=1ll*p*az%mod;}
(B[i-1]+=1ll*az*xs%mod)%=mod,(B[i-1]+=mod)%=mod;
}
}
int h[3010000];
int xa[10100],xc[10100];
int u[3001000];
void ot(int p,int c,int xs){
int mo=1;
for(int i=1;i<=c;i++)mo*=p;
memset(xa,0,sizeof(xa));memset(xc,0,sizeof(xc));
int t=d-1;
for(int i=0;i<d;i++)if(B[i]%p){t=i;break;}
int I=inv(B[t],mo),k=d-1-t;
xa[d-1]=I;
for(int i=0;i<d;i++)if(i!=k)xc[i]=(mo-1ll*B[d-1-i]*I%mo)%mo;
int L=d-1;
while(L>k&&!xc[L])L--;
while(L>k){
for(int i=L;i>k;i--)if(xc[i]){
(xa[i+d-1-k]+=1ll*xc[i]*I%mod)%=mo;
for(int j=0;j<d;j++)if(j!=k)(xc[i+(j-k)]+=1ll*(mo-1ll*B[d-1-j]*I%mo)%mo*xc[i]%mo)%=mo;
xc[i]=0;
}
L+=d;
while(L&&!xc[L])L--;
}
L=(d-1)*c;
while(L>=k&&!xa[L])L--;
I=inv(((1-xc[k])%mo+mo)%mo,mo);
for(int i=0;i<k;i++)xc[i]=1ll*xc[i]*I%mod;
for(int i=k;i<=L;i++)xa[i]=1ll*xa[i]*I%mod;
for(int i=0;i<m;i++){
int z=0;
for(int j=k;j<=L;j++)(z+=1ll*A[i+(j-k)]*xa[j]%mo)%=mo;
for(int j=0;j<k;j++)if(i>=k-j)(z+=1ll*xc[j]*u[i-(k-j)]%mo)%=mo;
u[i]=z,(h[i]+=1ll*u[i]*xs%mod)%=mod;
}
}
int main(){
scanf("%d%d%d%d",&n,&m,&d,&M);n*=d;
mod=M;
for(int i=2;i*i<=M;i++)if(!(M%i)){
t++,P[t]=i,C[t]=0;
while(!(M%i))M/=i,C[t]++;
}
if(M>1)t++,P[t]=M,C[t]=1;
for(int i=1;i<=t;i++){
int z=1;
for(int j=1;j<=C[i];j++)z*=P[i];
so(P[i],C[i],1ll*(mod/z)*inv(mod/z,z)%mod);
}
int pc=1;
for(int i=1;i<=t;i++)if(d%P[i])while(C[i]--)pc*=P[i];
int I=inv(B[0],pc),xs=1ll*(mod/pc)*inv(mod/pc,pc)%mod;
for(int i=0;i<m;i++){
h[i]=A[i];
for(int j=1;j<d;j++)if(j<=i)(h[i]-=1ll*B[j]*h[i-j]%pc)%=pc;
h[i]=1ll*h[i]*I%pc;
}
for(int i=0;i<m;i++)h[i]=(1ll*xs*h[i]%mod+mod)%mod;
for(int i=1;i<=t;i++)if(!(d%P[i])){
int z=1;for(int j=1;j<=C[i];j++)z*=P[i];
ot(P[i],C[i],1ll*(mod/z)*inv(mod/z,z)%mod);
}
int anss=0;
for(int i=0;i<m;i++)anss^=((h[i]+mod)%mod);
return printf("%d",anss),0;
}
