[SDOI2013]方程

Description
给定方程

\[X_1+X_2+...+X_n=M \]

我们对第\(1...n_1\)个变量进行一些限制:
\(X_1\leqslant A_1\)
\(X_2\leqslant A_2\)
\(X_{n_1}\leqslant A_{n_1}\)
我们对第\(n_1+1...n_1+n_2\)个变量进行一些限制:
\(X_{n_1+1}\geqslant A_{n_1+1}\)
\(X_{n_1+2}\geqslant A_{n_1+2}\)
\(X_{n_1+n_2}\geqslant A_{n_1+n_2}\)
求:在满足这些限制的前提下,该方程正整数解的个数。
答案可能很大,请输出对p取模后的答案,也即答案除以p的余数。

Input
输入含有多组数据,第一行两个正整数T,p。T表示这个测试点内的数据组数,p的含义见题目描述。
对于每组数据,第一行四个非负整数n,n1,n2,m。
第二行nl+n2个正整数,表示A1..n1+n2。请注意,如果n1+n2等于0,那么这一行会成为一个空行。

Output
共T行,每行一个正整数表示取模后的答案。

Sample Input
3 10007
3 1 1 6
3 3
3 0 0 5
3 1 1 3
3 3

Sample Output
3
6
0


如果不考虑前\(n_1\)个变量的限制,那么答案为

\[\binom{M-\sum\limits_{i=n_1+1}^{n_1+n_2}A_i-1}{n-1} \]

现在多了\(n_1\)个限制,但\(n_1\)很小,我们可以暴力容斥,枚举哪些变量不满足限制即可

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
#define Fi first
#define Se second
#define MK std::make_pair
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
typedef std::pair<int,int> pii;
inline char gc(){
	static char buf[1000000],*p1=buf,*p2=buf;
	return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
template<typename T>inline T frd(T x){
	int f=1; char ch=gc();
	for (;ch<'0'||ch>'9';ch=gc())	if (ch=='-')    f=-1;
	for (;ch>='0'&&ch<='9';ch=gc())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
template<typename T>inline T read(T x){
	int f=1;char ch=getchar();
	for (;ch<'0'||ch>'9';ch=getchar())	if (ch=='-')	f=-1;
	for (;ch>='0'&&ch<='9';ch=getchar())	x=(x<<1)+(x<<3)+ch-'0';
	return x*f;
}
inline void print(int x){
	if (x<0)    putchar('-'),x=-x;
	if (x>9)	print(x/10);
	putchar(x%10+'0');
}
template<typename T>inline T min(T x,T y){return x<y?x:y;}
template<typename T>inline T max(T x,T y){return x>y?x:y;}
template<typename T>inline T swap(T &x,T &y){T t=x; x=y,y=t;}
const int N=1e4;
namespace Math{
	int P[10],V[10],C[10],f[10][N+10],tot,SP;
	int mlt(int a,int b,int p){
		int res=1;
		for (;b;b>>=1,a=1ll*a*a%p)	if (b&1)	res=1ll*res*a%p;
		return res;
	}
	void prepare(int p){
		SP=p;
		for (int i=2;i*i<=p;i++){
			if (p%i)	continue;
			P[++tot]=i,V[tot]=1;
			while (p%i==0)	p/=i,V[tot]*=i,C[tot]++;
			f[tot][0]=1;
			for (int j=1;j<=V[tot];j++)	f[tot][j]=1ll*f[tot][j-1]*(j%i?j:1)%V[tot];
		}
		if (p!=1){
			f[++tot][0]=1;
			P[tot]=V[tot]=p,C[tot]=1;
			for (int j=1;j<=V[tot];j++)	f[tot][j]=1ll*f[tot][j-1]*(j%p?j:1)%V[tot];
		}
	}
	int gcd(int a,int b){return !b?a:gcd(b,a%b);}
	void exgcd(int a,int b,int &x,int &y){
		if (!b){x=1,y=0;return;}
		exgcd(b,a%b,x,y);
		int t=x; x=y,y=t-a/b*y;
	}
	int Ex_GCD(int a,int b,int c){
		int d=gcd(a,b),x,y;
		if (c%d)	return -1;
		a/=d,b/=d,c/=d;
		exgcd(a,b,x,y);
		x=(1ll*x*c%b+b)%b;
		return x;
	}
	pii work(int n,int i){
		if (n<=1)	return MK(1,0);
		int res=1ll*mlt(f[i][V[i]],n/V[i],V[i])*f[i][n%V[i]]%V[i];
		pii tmp=work(n/P[i],i); int cnt=n/P[i];
		return MK(1ll*res*tmp.Fi%V[i],tmp.Se+cnt);
	}
	int calc(int n,int m,int i){
		int res=1,cnt=0; pii tmp;
		tmp=work(n,i); cnt+=tmp.Se;
		res=1ll*res*tmp.Fi;
		
		tmp=work(m,i); cnt-=tmp.Se;
		res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
		
		tmp=work(n-m,i); cnt-=tmp.Se;
		res=1ll*res*Ex_GCD(tmp.Fi,V[i],1)%V[i];
		
		return cnt<C[i]?1ll*res*mlt(P[i],cnt,V[i])%V[i]:0;
	}
	int Ex_C(int n,int m){
		if (n<m)	return 0;
		if (tot==1)	return Ex_GCD(1,V[1],calc(n,m,1));
		int Ans=0;
		for (int i=1;i<=tot;i++)	Ans=(1ll*Ex_GCD(SP/V[i],V[i],1)*(SP/V[i])%SP*calc(n,m,i)%SP+Ans)%SP;
		return Ans;
	}
}
int main(){
	int T=read(0),p=read(0);
	Math::prepare(p);
	for (;T;T--){
		static int A[10];
		int n=read(0),x=read(0),y=read(0),m=read(0),Ans=0;
		for (int i=0;i<x;i++)	A[i]=read(0);
		for (int i=0;i<y;i++)	m-=read(0)-1;
		for (int sta=0;sta<1<<x;sta++){
			int del=0,cnt=0;
			for (int i=0;i<x;i++)	if (sta>>i&1)	cnt++,del+=A[i];
			Ans=(Ans+(cnt&1?-1:1)*Math::Ex_C(m-del-1,n-1))%p;
		}
		printf("%d\n",(Ans+p)%p);
	}
	return 0;
}
posted @ 2019-03-31 09:55  Wolfycz  阅读(349)  评论(0编辑  收藏  举报