【组合数取模】【SDOI2013】【BZOJ 3129】【JZOJ 3214】方程

Description

给定方程
$X_1+X_2+\dots +X_n=m$
我们对第$1\dots n_1$个变量 进行一些限制 ：
$X_1\le A_1$
$X_2\le A_2$
$\dots$
$X_{n1}\le A_{n1}$
我们对第$n_1+1,n_1+2\dots n_1+n_2$个变量进行一些限制 ：
$X_{n_1+1}\ge A_{n_1+1}$
$X_{n_1+2}\ge A_{n_1+2}$
$\dots$
$X_{n_1+n_2}\ge A_{n_1+n_2}$
求：在满足这些限制的前提下， 该方程正整数解的个数。
答案可能很大，请输出对 p取模 后的答案，也即答案除以p的余数。

Analysis

很好的一道题，综合了几个知识点，可以帮助复习算法（对于我这种蒟蒻只能学习了）。
有时要抽象题意，可是有时却要形象题意，这样的好处是更直观，例如此题，形象化：
有$m$个球，$n$个盒子，要把$m$个球放入$n$个盒子内，每个盒子至少放一个球，有些盒子还有必须放大于或小于某个球的限制，问方案数。
让我们一步一步来，从简单到复杂（不要说我生物学得太好）。
不考虑任何约束，易用隔板法证明，方案数为$C_{m-1}^{n-1}$。
对于那些$\ge A_i$的，我们可以让$m=m-(A_i-1)$，这个显然。
对于那些$\le A_i$的，很难搞，于是正难则反，反过来搞，变成$\ge$，发现要用容斥减重。
那么我们现在答案是一堆$C$的和。
是不是觉得能A了，激动了？
错！组合数涉及除法，模数不一定是素数。
但是，我们有
中国剩余定理！
中国剩余定理，又称孙子定理，是中国古人发现的神奇方法，被某人传到欧洲，因而得名。
这里我们的答案就是定理里面的$x$，将模数$p$分解成$m_1^{x_1}\ast m_2^{x_2}\cdots m_k^{x_k}$，对于每个$m_i$做一遍，求出对应的$a_i$，再合并。
所以，现在问题变成了如何计算$C_m^{n}modmo$
可以把$C_m^n$化成$a\ast m{o}^b$的形式，再$modmo$
为什么不可以直接模呢，因为若分母里有$mo$的倍数，那么互质条件不成立，逆元也就不能求了。
所以，我们的问题又化为更小的一个子问题，求$n!$里有多少个$mo$的及除了$mo$外的积。
直接暴力肯定会TLE，我们要找到快速方法。
设$fac[i]$表示前$i$个数不包括$mo$以外的乘积。这个只需要预处理到最大的$m{o}^b$就好了。
举个例子，$n=13,mo=3$，当前的ans为$dfs(n)$
$$n!=1\ast 2\ast 3\ast 4\ast 5\ast 6\ast 7\ast 8\ast 9\ast 10\ast 11\ast 12\ast 13$$
$$=1\ast 2\ast 4\ast 5\ast 7\ast 8\ast 10\ast 11\ast dfs(n/mo)\ast fac[nmodmo]$$
$$\equiv (1\ast 2{)}^4\ast dfs(n/mo)\ast fac[nmodmo]$$
递归处理即可。

Code

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
int n1,n2,tot,n,m,num,p,a[20];
ll fac[10210];
struct lyd
{
	int x,y,z;
}c[10];
ll qmi(int x,int n,int mo)
{
	ll t=1;
	for(;n;n>>=1)
	{
		if(n&1) t=t*x%mo;
		x=x*x%mo;
	}
	return t;
}
void exgcd(int a,int b,ll &x,ll &y)
{
	if(!b)
	{
		x=1,y=0;
		return;
	}
	ll x1,y1;
	exgcd(b,a%b,x1,y1);
	x=y1,y=x1-a/b*y1;
}
ll ny(int n,int mo)
{
	ll x,y;
	exgcd(n,mo,x,y);
	return (x%mo+mo)%mo;
}
void pri(int n)
{
	num=0;
	for(int i=2;i*i<=n;i++)
		if(n%i==0)
		{
			c[++num].x=i,c[num].y=0;
			while(n%i==0) n/=i,c[num].y++;
		}
	if(n!=1) c[++num].x=n,c[num].y=1;
}
ll calcfac(int n,int k,int hjy)
{
	if(n<c[k].x) return fac[n];
	int mo=c[k].z;
	tot+=hjy*(n/c[k].x);
	return calcfac(n/c[k].x,k,hjy)*qmi(fac[mo],n/mo,mo)%mo*fac[n%mo]%mo;
}
ll CC(int m,int n,int k)
{
	tot=0;
	ll t=1;
	int mo=c[k].z;
	t=calcfac(m,k,1)*ny(calcfac(n,k,-1)*calcfac(m-n,k,-1)%mo,mo)%mo;
	t=t*qmi(c[k].x,tot,mo)%mo;
	return t;
}
ll C(int m,int n)
{
	if(m<n) return 0;
	pri(p);
	fo(i,1,num) c[i].z=qmi(c[i].x,c[i].y,p+1);
	ll ans=0;
	fo(i,1,num)
	{
		fac[0]=1;
		fo(j,1,c[i].z)
		{
			fac[j]=fac[j-1];
			if(j%c[i].x!=0) fac[j]=fac[j]*j%c[i].z;
		}
		(ans+=(p/c[i].z)*ny(p/c[i].z,c[i].z)*CC(m,n,i))%=p;
	}
	return ans;
}
int main()
{
	int _;
	ll m1;
	scanf("%d %d",&_,&p);
	while(_--)
	{
		scanf("%d %d %d %d",&n,&n1,&n2,&m1);
		fo(i,1,n1+n2) scanf("%d",&a[i]);
		fo(i,n1+1,n1+n2) m1-=a[i]-1;
		ll ans=C(m1-1,n-1);
		fo(i,1,(1<<n1)-1)
		{
			ll m=m1,t=0;
			fo(j,0,n1-1)
				if(i&(1<<j)) t++,m-=a[j+1];
			t=(t&1)?(-1):(1);
			ans=((ans+t*C(m-1,n-1)%p)%p+p)%p;
		}
		printf("%lld\n",ans);
	}
}