BZOJ 3129 方程

发表于 | 更新于

Link

就是这么一道题让我搭上了几乎一天。。

Solution

个人感觉两步都很难啊。 首先，给定∑xi=N\sum x_i=N这样的方程，如果是求正整数解的方案数目是很好求的。 对于xi≥aix_i\ge a_i的限制，可以直接把总和−(ai−1)-(a_i-1)，相当于把ai−1a_i-1钦定给了xix_i，就可以转化为正整数解了。 对于xi≤aix_i\le a_i的限制就诡了。正解是容斥，设f(...)f(...)为钦定......不满足条件的方案数目，那么Ans=f()−∑f(i)+∑∑f(i,j)−∑∑∑f(i,j,k)...Ans=f()-\sum f(i)+\sum\sum f(i,j)-\sum\sum\sum f(i,j,k)... 之前没怎么研究过容斥的题目，感觉这个条件给的很玄。 于是我查了查资料，翻了翻书，在《组合数学》里的一句话启发了我：

AiA_i是SS中具有性质PiP_i(也可能还具有其他的一些性质)的子集

这样，某个xix_i一定不对就是一个这样的性质。 至于第二步，就是这篇文章了；

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=10;
typedef long long LL;
int P;
inline int Pow(int a,int v)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res*=a;
		a*=a;
		v>>=1;
	}
	return res;
}
inline int Pow(int a,int v,int P)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res=(LL)res*a%P;
		a=(LL)a*a%P;
		v>>=1;
	}
	return res;
}
int exgcd(int a,int b,LL &x,LL &y)
{
	if(!b)
	{
		x=1,y=0;
		return a;
	}
	else
	{
		int d=exgcd(b,a%b,x,y);
		LL nx=y,ny=x-(a/b)*y;
		x=nx,y=ny;
		return d;
	}
}
inline int inv(int a,int n)
{
	LL x,y;exgcd(a,n,x,y);
	return (x%n+n)%n;
}
namespace iLucas
{
	int P,dc;
	int dv[MAXN],dt[MAXN],dvt[MAXN];
	//A[]
	//M[i]=M/m[i]
	//t[i]*M[i]=1 (mod m[i])
	int A[MAXN],M[MAXN],t[MAXN],m[MAXN],tM[MAXN];
	const int FN=1e5,MAXF=FN+10;
	int table[MAXN][MAXF];
	inline void Add(int np,int &x)
	{
		++dc;dv[dc]=np;dt[dc]=0,dvt[dc]=1;
		while(x%np==0) dt[dc]++,x/=np,dvt[dc]*=np;
		int *ptr=table[dc];ptr[0]=1;
		for(int i=1;i<=dvt[dc];i++)
			if(i%np==0) ptr[i]=ptr[i-1];
			else ptr[i]=(LL)ptr[i-1]*i%dvt[dc];
		m[dc]=dvt[dc];M[dc]=P/m[dc];t[dc]=inv(M[dc],m[dc]);
		tM[dc]=(LL)t[dc]*M[dc]%P;
	}
	inline void Init(int x)
	{
		P=x;
		for(int i=2;i*i<=x;i++)
			if(x%i==0) Add(i,x);
		if(x!=1) Add(x,x);
	}
	int Calc(int n,int idx)
	{
		int *ptr=table[idx],v=dv[idx],vt=dvt[idx];
		if(n<v) return ptr[n];
		return (LL)Calc(n/v,idx)*Pow(ptr[vt],n/vt,vt)*ptr[n%vt]%vt;
	}
	int C(int n,int m,int idx)
	{
		int times,v,t,vt;
		v=dv[idx],t=dt[idx],vt=dvt[idx];times=0;
		for(int i=n;i;i/=v) times+=i/v;
		for(int i=m;i;i/=v) times-=i/v;
		for(int i=n-m;i;i/=v) times-=i/v;
		if(times>=t) return 0;
		else
		{
			int pt,fn,fm,fnmm;
			pt=Pow(v,times);
			fn=Calc(n,idx),fm=Calc(m,idx),fnmm=Calc(n-m,idx);
			return (LL)pt*fn*inv(fm,vt)*inv(fnmm,vt)%vt;
		}
	}
	inline int CRT()
	{
		int res=0;
		for(int i=1;i<=dc;i++) res=((LL)A[i]*tM[i]+res)%P;
		return res;
	}
	inline int C(int n,int m)
	{
		if(n<m) return 0;
		for(int i=1;i<=dc;i++) A[i]=C(n,m,i);
		return CRT();
	}
}
int a[MAXN],ac,n,m;
//前ac个数<=a[i] 所有的数>=1的方案数
void WORK()
{
	int n1,n2,t;
	get(n),get(n1),get(n2),get(m);
	for(int i=1;i<=n1;i++) get(a[i]);
	ac=n1;
	for(int i=1;i<=n2;i++) get(t),m-=t-1;
	int sum,cnt,ANS=0,res;
	for(int i=(1<<n1)-1;~i;i--)
	{
		sum=m,cnt=0;
		for(int j=1;j<=ac;j++)
			if(i&(1<<j>>1)) cnt++,sum-=a[j];
		res=iLucas::C(sum-1,n-1);
		if(cnt&1) ANS-=res,((ANS%=P)+=P)%=P;
		else ANS+=res,ANS%=P;
	}
	printf("%d\n",ANS);
}
int main()
{
	int T;get(T),get(P);iLucas::Init(P);
	while(T--) WORK();
	return 0;
}
/* AC Record(Bugs) *
 * WA baoint +=
 * WA invine 
 * * * * * * * * * */