[loj3075]组合数求和

Subtask1：​​​$m,ndle 2 imes 10^{3}$

对$M$质因数分解，假设$M=prod_{i=1}^{k}p_{i}^{alpha_{i}}$（其中$p_{i}$为素数），对每个$i$求出$f(j) mod p_{i}^{alpha_{i}}$的值，再通过exgcd即可求出$f(j)$，注意到$sum log p_{i}^{alpha_{i}}=log M$，即exgcd的复杂度仅为$o(mlog M)$

在模$p^{alpha}$的意义下，只需要将数值以$rcdot p^{k}$的形式描述（其中$p otmid r$），即可支持求逆元

预处理$nd$以内阶乘即逆元，复杂度为$o(ndlog M)$，即可$o(1)$求出${nchoose m}mod p^{alpha}$，那么暴力计算即可

时间复杂度为$o(nmlog M)$，可以通过

Subtask2：$d=1$

对于多项式$F(x)$，记$F(x)[x^{n}]$表示该多项式$n$次项系数

根据二项式定理，有${nchoose m}=(1+x)^{n}[x^{m}]$，即$f(j)=(sum_{i=0}^{n-1}(1+x)^{id})[x^{j}]mod M$

记后者为$C(x)$，根据等比数列求和，即$C(x)=frac{(1+x)^{nd}-1}{(1+x)^{d}-1}mod x^{m}$

由于分子和分母同时除以$x$，即分别为$egin{cases}A(x)=sum_{i=1}^{nd}{ndchoose i}x^{i-1}\B(x)=sum_{i=1}^{d}{dchoose i}x^{i-1}end{cases}$，那么$C(x)=frac{A(x)}{B(x)}mod x^{m}$

由于$d=1$，即$B(x)=1$，那么$C(x)=A(x) mod x^{m}$

考虑$A(x)[x^{i}]={ndchoose i+1}$，将$i$从小到大枚举，每次即乘上$frac{nd-i}{i+1}$，关于如何支持除法与Subtask1相同（对$M$质因数分解并将数值以$rcdot p^{k}$的形式描述），即可$o(mlog M)$求出$A(x)$

时间复杂度为$o(mlog M)$，可以通过

Subtask3：$mle 8 imes 10^{3}$，$M=998244353$

当$M=998244353$时，$A(x)$可以线性预处理逆元做到$o(m)$，$B(x)$多项式求逆即可

时间复杂度为$o(m log m)$，可以通过

（由于本做法与正解关系不大，这里就不实现了）

Subtask4：$gcd(d,M)=1$

仍考虑$C(x)=frac{A(x)}{B(x)}mod x^{m}$，也即$B(x)cdot C(x)equiv A(x)(mod x^{m})$

考虑$i$次项系数，即$forall 0le i<m,sum_{j=0}^{i}B(x)[x^{j}]cdot C(x)[x^{i-j}]=A(x)[x^{i}]$

再代入具体的式子，即$sum_{j=0}^{i}{dchoose j+1}cdot C(x)[x^{i-j}]=A(x)[x^{i}]$

将$j e 0$的部分减到右边并同除以$d$，即$C(x)[x^{i}]=frac{A(x)[x^{i}]-sum_{j=1}^{i}{dchoose j+1}cdot C(x)[x^{i-j}]}{d}$

$A(x)$可以$o(mlog M)$求出，而求和式在$jge d$时无意义，可以$o(d)$求出

由于$gcd(d,M)=1$，也即$p otmid d$，可以扩欧求出$d$在模$p^{alpha}$意义下的逆元来支持除法

时间复杂度为$o(m(d+log M))$，可以通过

Subtask5：$d$为质数

当$p otmid d$时，仍用逆元的方式处理即可，以下考虑$pmid d$的情况：

由于$d$是素数，那么$pmid d$当且仅当$p=d$

此时，将原式变形，可得$C(x)[x^{i}]=A(x)[x^{i+d-1}]-sum_{j=1}^{d-1}{dchoose j}C(x)[x^{i+j}]$

（虽然这个等式仅在$i+dle m$时成立，但令其在$i$更大时成立不影响正确性）

不断迭代，即按照$i$从小到大维护当前$C(x)[x^{i}]$的系数，并将其乘上$-{dchoose j}$转移到$C(x)[x^{i+j}]$

每转移一次都会让$p$的幂次增加1，那么最多转移$alpha$次即会在模$p^{alpha}$意义下为0

由此，我们发现$C(x)[x^{i}]$可以被描述为关于$A(x)[x^{i+d-1}],A(x)[x^{i+d}],...,A(x)[x^{i+dalpha}]$的一个线性的式子，在一开始迭代$o(d^{2}alpha)$预处理出来，再利用此递推式$o(mdalpha)$计算即可

时间复杂度为$o(mdlog M)$，可以通过

Subtask6：无特殊限制

仍考虑$p otmid d$的情况，此时前面的做法并不一定正确，原因是存在$1le j<d$使得${dchoose j}$使得不是$p$的倍数

对于这个问题，令$t$为最大的$j$使得${dchoose j}$不是$p$的倍数（其中$1le j<d$），将原式变形，即
$$
C(x)[x^{i+t}]=frac{A(x)[x^{i+d-1}]-sum_{j=0}^{t-1}{dchoose j}C(x)[x^{i+j}]-sum_{j=t+1}^{d-1}{dchoose j}C(x)[x^{i+j}]}{{dchoose t}}
$$
注意到每向后迭代一次（最后一项），都会让$p$的幂次增加1

换言之，我们只需要将其分层迭代，每一层按次数从大到小迭代，并且当迭代到$i$之前后就不需要再迭代（直接取该值即可），那么每一层要迭代$o(dalpha)$个数，每一个数迭代复杂度为$o(d)$，最多$alpha$层，复杂度即$o(d^{2}alpha^{2})$

同样预处理出这个线性的式子，再$o(mdalpha)$计算即可

时间复杂度为$O(d^{2}log^{2}M+mdlog M)$，可以通过

（事实上，前面Subtask4和Subtask5即$t=d-1$和$t=0$的情况）

#include<bits/stdc++.h>
using namespace std;
#define N 3000005
#define D 105
#define L 35
vector<pair<int,int> >v;
int n,m,d,M,p,s,mod,Ans,mi[L],C[D][D],A[N+D*L],tot[D*L],Tot[D*L],f[N],ans[N];
int exgcd(int a,int b,int &x,int &y){
    if (!b){
        x=1,y=0;
        return a;
    }
    int g=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return g;
}
int get_inv(int k){
    int x,y;
    exgcd(k,mod,x,y);
    return (x%mod+mod)%mod;
}
struct num{
    int r,k;
    num(){
        r=1,k=0;
    }
    num(int rr){
        r=rr,k=0;
        if (!r)return;
        while (r%p==0){
            r/=p;
            k++;
        }
    }
    num(int rr,int kk){
        r=rr,k=kk;
    }
    num operator * (const num &a)const{
        return num(1LL*r*a.r%mod,k+a.k);
    }
    num inv(){
        return num(get_inv(r),-k);
    }
    int get(){
        if (k>=s)return 0;
        return 1LL*r*mi[k]%mod;
    }
};
void calc(){
    mi[0]=1;
    for(int i=1;i<s;i++)mi[i]=mi[i-1]*p;
    for(int i=0;i<=d;i++){
        C[i][0]=C[i][i]=1;
        for(int j=1;j<d;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
    num ans=num(1);
    for(int i=0;i<=m+s*d;i++){
        if (n*d<i)ans=num(0);
        else ans=ans*num(n*d-i)*num(i+1).inv();
        A[i]=ans.get();
    }
    int t=0;
    for(int i=1;i<d;i++)
        if (C[d][i]%p)t=i;
    int inv=get_inv(C[d][t]);
    memset(tot,0,sizeof(tot));
    memset(Tot,0,sizeof(Tot));
    tot[d]=1;
    for(int i=0;i<s;i++)
        for(int j=d*s;j>=d;j--){
            Tot[j]=(Tot[j]+1LL*inv*tot[j])%mod;
            for(int k=0;k<d;k++)
                if (k!=t)tot[j+k-t]=(tot[j+k-t]-1LL*C[d][k]*inv%mod*tot[j]%mod+mod)%mod;
            tot[j]=0;
        }
    for(int i=0;i<m;i++){
        f[i]=0;
        for(int j=1;j<=min(i,d-1);j++)f[i]=(f[i]+1LL*tot[d-j]*f[i-j]%mod)%mod;
        for(int j=d;j<=d*s;j++)f[i]=(f[i]+1LL*Tot[j]*A[i+j-t-1])%mod;
    }
}
int main(){
    scanf("%d%d%d%d",&n,&m,&d,&M);
    int k=M;
    for(int i=2;i*i<=k;i++){
        int s=0;
        while (k%i==0){
            k/=i;
            s++;
        }
        if (s)v.push_back(make_pair(i,s));
    }
    if (k>1)v.push_back(make_pair(k,1));
    k=1;
    for(int i=0;i<v.size();i++){
        p=v[i].first,s=v[i].second,mod=1;
        for(int j=0;j<s;j++)mod*=p;
        calc();
        int x=get_inv(k),kk=k*mod;
        for(int j=0;j<m;j++)ans[j]=(1LL*x*(f[j]-ans[j]+kk)%kk*k+ans[j])%kk;
        k=kk;
    }
    for(int i=0;i<m;i++)Ans^=ans[i];
    printf("%d",Ans); 
}