NKOJ 2640 （SDOI 2013）方程（扩展Lucas+容斥原理）

发表于 |
字数统计: 943字 | 阅读时长 ≈ 5min

P2640【SDOI2013 R1 Day2】方程

问题描述

这里写图片描述

输入格式

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

输出格式

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

样例输入

3 10007
3 1 1 6
3 3
3 0 0 5
3 1 1 3
3 3

样例输出

3
6
0

提示

这里写图片描述

考虑没有限制的情况，那么显然用隔板法求解，$Ans=C_{m-1}^{n-1}$
考虑大于等于的限制，那么显然直接$m=m-\sum A_i-1$即可。
考虑小于等于的限制，不好做，但限制数很少，考虑容斥。
那么$Ans=C_{m-sum}^{n-1}-有一个X_i不满足条件的方案+有两个X_i不满足条件的方案……$
显然，后面的东西也是一个隔板法，那么剩下的问题就是算组合数。

如果p是质数，直接Lucas，而本题p可能是合数，需要用到扩展Lucas，大体思路是利用CRT解方程，
即将p分解质因数，分别对每个$p_i^{a_i}$计算答案得到一组同余方程，然后因为模数互质，可以用CRT合并得到最终答案，具体请百度扩展Lucas

代码：

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
#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
#define N 15
#define ll long long
using namespace std;
struct node{ll a,b;};
bool operator<(node a,node b)
{
	if(a.b==b.b)return a.a<b.a;
	return a.b<b.b;
}
map<node,ll>Q;
ll T,p,A[N];
ll QM(ll a,ll b,ll c)
{
	ll o=1;
	while(b)
	{
		if(b&1)o=o*a%c;
		b>>=1;a=a*a%c;
	}
	return o;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)x=1,y=0;
	else exgcd(b,a%b,y,x),y-=a/b*x;
}
ll inv(ll a,ll b)
{
	ll x,y;exgcd(a,b,x,y);
	x=(x%b+b)%b;
	return x?x:x+b;
}
ll fac(ll n,ll px,ll pi)
{
	if(!n)return 1;
	node tmp=(node){n,px};
	if(Q[tmp])return Q[tmp];
	ll i,ans=1;
	if(n/px)
	{
		for(i=2;i<=px;i++)if(i%pi)ans=ans*i%px;
		ans=QM(ans,n/px,px);
	}
	for(i=2;i<=n%px;i++)if(i%pi)ans=ans*i%px;
	return Q[tmp]=ans*fac(n/pi,px,pi)%px;
}
ll C(ll n,ll m,ll px,ll pi)
{
	if(m>n)return 0ll;
	ll a=fac(n,px,pi),b=fac(m,px,pi),c=fac(n-m,px,pi);
	ll i,k=0,ans=0;
	for(i=n;i;i/=pi)k+=i/pi;
	for(i=m;i;i/=pi)k-=i/pi;
	for(i=n-m;i;i/=pi)k-=i/pi;
	ans=a*inv(b,px)%px*inv(c,px)%px*QM(pi,k,px)%px;
	return ans*(p/px)%p*inv(p/px,px)%p;
}
ll Cal(ll n,ll m)
{
	ll i,ans=0,x=p,px;
	for(i=2;i<=x;i++)
	{
		if(x%i)continue;px=1;
		while(x%i==0)px*=i,x/=i;
		ans=(ans+C(n,m,px,i))%p;
	}
	return ans;
}
int main()
{
	ll n,n1,n2,m,i,j,k,x,y,ans,S;
	scanf("%lld%lld",&T,&p);
	while(T--)
	{
		scanf("%lld%lld%lld%lld",&n,&n1,&n2,&m);
		for(i=1;i<=n1;i++)scanf("%lld",&A[i]);
		for(i=1;i<=n2;i++)scanf("%lld",&x),m-=x-1;
		ans=Cal(m-1,n-1);S=(1<<n1)-1;
		for(i=1;i<=S;i++)
		{
			x=0;y=m;
			for(j=1;j<=n1;j++)if(i>>j-1&1)x++,y-=A[j];
			if(x&1)ans-=Cal(y-1,n-1),ans%=p;
			else ans+=Cal(y-1,n-1),ans%=p;
		}
		printf("%lld\n",(ans%p+p)%p);
	}
}