BZOJ4338 BJOI2015糖果

题目描述

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\le N,M\le 10^5,1\le P,K\le 2\times 10^9$

题解

每一列单调不减，即为在$k$个数中选$m$个数的可重组合数，即$C_{k+m-1}^k$，令$r=C_{k+m-1}^k$。

有$n$行，每行互不相同，那么总共有$A_r^n$种方案，那么答案就是$A_r^n\%p$。

考虑如何求这个数。我们发现$A_r^n\%p=(A_{r\%p}^{n\%p})\%p$，所以我们可以先求$r\%p$，即$C_{k+m-1}^k\%p$。

假设我们要求$C_n^m\%p$。设$p=p_1^{r_1}{p}_2^{r_2}\cdots p_k^{r_k}$，其中$p_i$为质数。我们可以先求出$C_n^m\%p_1^{r_1},C_n^m\%p_2^{r_2},\dots ,C_n^m\%p_k^{r_k}$的值$a_1,a_2,\dots ,a_k$。

我们把$C_n^m$看作未知数$x$，可以得到以下方程组：

$$\{\begin{array}{c}x\equiv a_1(\mathrm{mod}{p}_1^{r_1})\\ x\equiv a_2(\mathrm{mod}{p}_2^{r_2})\\ x\equiv a_3(\mathrm{mod}{p}_3^{r_3})\\ ......\\ x\equiv a_n(\mathrm{mod}{p}_k^{r_k})\end{array}$$

利用中国剩余定理，我们可以求出$x$，它是以$p$为周期出现的无穷多个解。而在$[0,p)$这个周期的解，就是$C_n^m\%p$后的值。

那么$a_1,a_2\dots ,a_k$怎么求呢？

$a_1=C_n^m\%p_1^{r_1}=\frac{n^{\underline{m}}}{m!}\%p_1^{r_1}$

对于上面这个式子，我们可以将分子和分母的质因子$p_1$个数求出来。因为$C_n^m$是一个整数，所以分子含有的$p_1$的数量一定大于等于分母含有的$p_1$的数量。将分母中的$p_1$约分去掉，此时的分母不含质因数$p_1$，也就是与$p_1^{r_1}$互质，那么就可以用欧拉定理求分母的逆元，这样即可求出$a_1$的值。

求出$C_{k+m-1}^k\%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;
}