BZOJ3129 [Sdoi2013]方程 【扩展Lucas】

题目

给定方程
X1+X2+. +Xn=M
我们对第l..N1个变量进行一些限制：
Xl < = A
X2 < = A2
Xn1 < = An1
我们对第n1 + 1..n1+n2个变量进行一些限制：
Xn1+l > = An1+1
Xn1+2 > = An1+2
Xnl+n2 > = Anl+n2
求：在满足这些限制的前提下，该方程正整数解的个数。
答案可能很大，请输出对p取模后的答案，也即答案除以p的余数。

输入格式

输入含有多组数据，第一行两个正整数T，p。T表示这个测试点内的数据组数，p的含义见题目描述。
对于每组数据，第一行四个非负整数n，n1，n2，m。
第二行nl+n2个正整数，表示A1..n1+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

【样例说明】

对于第一组数据，三组解为(1，3，2)，(1,4,1),(2,3,1)

对于第二组数据，六组解为(1，1，3)，(1,2,2),(1,3,1),(2,1,2),(2,2,1),(3,1,1)

提示

n < = 10^9 , n1 < = 8 , n2 < = 8 , m < = 10^9 ,p<=437367875

对于l00%的测试数据： T < = 5，1 < = A1..n1_n2 < = m，n1+n2 < = n

题解

组合数经典套路：
对于下限的限制，我们预先分配那么多的数，就去掉了这个下限
对于上限的限制，我们容斥一下哪些数超过了限制，就去掉了上限

被卡常卡哭了，，
膜了一下别人的代码才发现扩展Lucas在计算前预处理一下阶乘会快很多

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
#define res register
using namespace std;
const int maxn = 11000,maxm = 100005,INF = 1000000000;
inline int read(){
	int out = 0,flag = 1; char c = getchar();
	while (c < 48 || c > 57){if (c == '-') flag = -1; c = getchar();}
	while (c >= 48 && c <= 57){out = (out << 3) + (out << 1) + c - 48; c = getchar();}
	return out * flag;
}
int md,A[11];
int pi[11],pk[11],cnt,fac[maxn];
void pre(int pi,int pk){
	fac[0] = 1;
	for (int i = 1; i <= pk; i++)
		if (i % pi) fac[i] = 1ll * fac[i - 1] * i % pk;
		else fac[i] = fac[i - 1];
}
void init(){
	int x = md;
	for (res int i = 2; i * i <= x; i++)
		if (x % i == 0){
			++cnt; pi[cnt] = i; pk[cnt] = 1;
			while (x % i == 0) pk[cnt] *= i,x /= i;
		}
	if (x - 1) ++cnt,pi[cnt] = pk[cnt] = x;
}
inline int qpow(int a,int b,int md){
	int ans = 1;
	for (; b; b >>= 1,a = 1ll * a * a % md)
		if (b & 1) ans = 1ll * ans * a % md;
	return ans % md;
}
void exgcd(LL a,LL b,LL& d,LL& x,LL &y){
	if (!b) {x = 1; y = 0; d = a;}
	else exgcd(b,a % b,d,y,x),y -= a / b * x;
}
inline LL inv(LL a,LL P){
	if (!a) return 0;
	LL d,x,y; exgcd(a,P,d,x,y);
	x = (x % P + P) % P; if (!x) x += P;
	return x;
}
inline LL Fac(LL n,LL pi,LL pk){
	if (!n) return 1;
	int ans = qpow(fac[pk],n / pk,pk);
	return 1ll * ans * fac[n % pk] % pk * Fac(n / pi,pi,pk) % pk;
}
inline int C(int n,int m,int pi,int pk){
	if (m > n) return 0;
	int a = Fac(n,pi,pk),b = Fac(m,pi,pk),c = Fac(n - m,pi,pk),ans,k = 0;
	for (res int i = n; i; i /= pi) k  += i / pi;
	for (res int i = m; i; i /= pi) k  -= i / pi;
	for (res int i = n - m; i; i /= pi) k  -= i / pi;
	ans = 1ll * a * inv(b,pk) % pk * inv(c,pk) % pk * qpow(pi,k,pk) % pk;
	return 1ll * ans * (md / pk) % md * inv(md / pk,pk) % md;
}
inline int exlucas(int n,int m){
	if (m > n) return 0;
	int ans = 0;
	for (res int i = 1; i <= cnt; i++)
		pre(pi[i],pk[i]),ans = (ans + C(n,m,pi[i],pk[i])) % md;
	return ans;
}
void solve(){
	int N = read(),n1 = read(),n2 = read(),M = read();
	M -= N;
	for (res int i = 1; i <= n1; i++) A[i] = read();
	for (res int i = 1; i <= n2; i++) M -= read() - 1;
	int maxv = (1 << n1) - 1;
	int re = 0;
	for (res int s = 0; s <= maxv; s++){
		int m = M,pos = 1;
		for (res int i = 1,t = s; i <= n1; i++, t >>= 1)
			if (t & 1) m -= A[i],pos = -pos;
		if (m + N - 1 < 0) continue;
		re = (re + 1ll * pos * exlucas(m + N - 1,N - 1) % md) % md;
	}
	re = (re % md + md) % md;
	printf("%d\n",re);
}
int main(){
	int T = read(); md = read();
	init();
	while (T--) solve();
	return 0;
}