题目描述
Alice 正在教她的弟弟 Bob 学数学。
每天,Alice 画一个N行M 列的表格,要求 Bob在格子里填数。
Bob已经学会了自然数1到K的写法。因此他在每个格子里填1 ~ K之间的整数。
Alice 告诉 Bob,如果 Bob 填写完表格的 N*M 个数以后,每行的数从第 1 列到第 M
列单调不减,并且任意两行至少有一列的数不同,而且以前 Bob 没有填写过相同的表格,
那么Alice 就给Bob吃一颗糖果。
Bob想知道,如果每天填写一遍表格,最多能吃到多少颗糖果。
答案模P输出。
输入
第一行,四个整数依次是N, M, K, P。
输出
输出一行,一个整数,表示答案模P 后的结果。
样例
样例输入1
1 3 3 10
样例输出1
0
样例输入2
2 2 2 10
样例输出2
6
样例解释
样例1:表格只有一行。每个格子可以填写1 ~ 3。有10种填写方法,依次为 1 1 1, 1 1 2, 1 1 3, 1 2 2, 1 2 3, 1 3 3, 2 2 2, 2 2 3, 2 3 3, 3 3 3。
样例2: 表格有两行,有6种填写方法,依次为 1 1/1 2, 1 1/2 2, 1 2/1 1, 1 2/2 2, 2 2/1 1, 2 2/1 2。
数据范围
1≤N,M≤105,1≤P,K≤2×1091\leq N,M\leq 10^5,1\leq P,K\leq 2\times 10^91≤N,M≤105,1≤P,K≤2×109
题解
每一列单调不减,即为在kkk个数中选mmm个数的可重组合数,即Ck+m−1kC_{k+m-1}^kCk+m−1k,令r=Ck+m−1kr=C_{k+m-1}^kr=Ck+m−1k。
有nnn行,每行互不相同,那么总共有ArnA_r^nArn种方案,那么答案就是Arn%pA_r^n\% pArn%p。
考虑如何求这个数。我们发现Arn%p=(Ar%pn%p)%pA_r^n\% p=(A_{r\%p}^{n\% p})\% pArn%p=(Ar%pn%p)%p,所以我们可以先求r%pr\%pr%p,即Ck+m−1k%pC_{k+m-1}^k\%pCk+m−1k%p。
假设我们要求Cnm%pC_n^m\% pCnm%p。设p=p1r1p2r2⋯pkrkp=p_1^{r_1}p_2^{r_2}\cdots p_k^{r_k}p=p1r1p2r2⋯pkrk,其中pip_ipi为质数。我们可以先求出Cnm%p1r1,Cnm%p2r2,…,Cnm%pkrkC_n^m\%p_1^{r_1},C_n^m\%p_2^{r_2},\dots,C_n^m\%p_k^{r_k}Cnm%p1r1,Cnm%p2r2,…,Cnm%pkrk的值a1,a2,…,aka_1,a_2,\dots,a_ka1,a2,…,ak。
我们把CnmC_n^mCnm看作未知数xxx,可以得到以下方程组:
{x≡a1(modp1r1)x≡a2(modp2r2)x≡a3(modp3r3)......x≡an(modpkrk) \left\{ \begin{matrix} x\equiv a_1\pmod{p_1^{r_1}}\\ x\equiv a_2\pmod{p_2^{r_2}}\\ x\equiv a_3\pmod{p_3^{r_3}}\\ ......\\ x\equiv a_n\pmod{p_k^{r_k}} \end{matrix} \right. ⎩⎨⎧x≡a1(modp1r1)x≡a2(modp2r2)x≡a3(modp3r3)......x≡an(modpkrk)
利用中国剩余定理,我们可以求出xxx,它是以ppp为周期出现的无穷多个解。而在[0,p)[0,p)[0,p)这个周期的解,就是Cnm%pC_n^m\%pCnm%p后的值。
那么a1,a2…,aka_1,a_2\dots,a_ka1,a2…,ak怎么求呢?
a1=Cnm%p1r1=nm‾m!%p1r1a_1=C_n^m\%p_1^{r_1}=\dfrac{n^{\underline{m}}}{m!}\%p_1^{r_1}a1=Cnm%p1r1=m!nm%p1r1
对于上面这个式子,我们可以将分子和分母的质因子p1p_1p1个数求出来。因为CnmC_n^mCnm是一个整数,所以分子含有的p1p_1p1的数量一定大于等于分母含有的p1p_1p1的数量。将分母中的p1p_1p1约分去掉,此时的分母不含质因数p1p_1p1,也就是与p1r1p_1^{r_1}p1r1互质,那么就可以用欧拉定理求分母的逆元,这样即可求出a1a_1a1的值。
求出Ck+m−1k%pC_{k+m-1}^k\%pCk+m−1k%p,然后即可求出答案。
code
#include<bits/stdc++.h>
using namespace std;
int tot=0,p[105],q[105];
long long n,m,k,mod,ny,jc,vt=1,x,y,now=0,ans=1,a[105],r[105];
long long mi(long long t,long long v){
if(v==0) return 1;
long long re=mi(t,v/2);
re=re*re%mod;
if(v&1) re=re*t%mod;
return re;
}
void exgcd(long long c,long long d){
if(d==0){
x=1;y=0;
return;
}
exgcd(d,c%d);
long long t=x;x=y;y=t-c/d*y;
}
int main()
{
scanf("%lld%lld%lld%lld",&n,&m,&k,&mod);
long long v=mod;
for(int i=2;i*i<=v;i++){
if(v%i==0){
p[++tot]=i;r[tot]=1;
while(v%i==0){
++q[tot];r[tot]*=i;
v/=i;
}
}
}
if(v>1){
p[++tot]=r[tot]=v;q[tot]=1;
}
for(int i=1;i<=tot;i++){
long long v=0;
ny=jc=1;
for(int j=1;j<=m;j++){
long long t=j;
while(t%p[i]==0) t/=p[i],--v;
ny=ny*t%r[i];
t=(m+k-1)-j+1;
while(t%p[i]==0) t/=p[i],++v;
jc=jc*t%r[i];
}
a[i]=mi(ny,r[i]-r[i]/p[i]-1)%r[i]*jc%r[i]*(mi(p[i],v)%r[i])%r[i];
vt=vt*r[i];
}
for(int i=1;i<=tot;i++){
exgcd(vt/r[i],r[i]);
x=(x%r[i]+r[i])%r[i];
now=(now+vt/r[i]*a[i]*x%vt)%vt;
}
n%=mod;
for(long long i=now;i>=now-n+1;i--) ans=ans*i%mod;
printf("%lld",ans);
return 0;
}