BZOJ3129/洛谷P3301方程(SDOI2013)容斥原理+扩展Lucas定理

时间:2021-06-01 19:55:11

题意:给定方程x1+x2+....xn=m,每个x是正整数。但是对前n1个数做了限制x1<=a1,x2<=a2...xn1<=an1,同时对第n1+1到n1+n2个数也做了限制xn1+1>=an1+1....xn1+n2>=an1+n2,输出方程解个数。

解法:首先如果对数字没有任何要求(应该是只要求是非负数)的话,答案就是C(n+m+1,m+1)原理是隔板法。但是此题有各种限制,我们想办法解决限制使得答案往无限制上面靠。

首先是解决要正整数,那么每个数字减一即可,就是m-=n。

然后对于n1+1到n1+n2的限制,大于等于的限制也简单,也是xi>=ai的话就使xi-=ai即可。

好了到最麻烦的xi<=ai限制,这个限制我们没办法通过数字处理来解决。但是我们观察到这些限制个数小于等于8个,自然会想到容斥原理。

那么容斥原理常规套路,无限制-至少一个突破限制(也就是xi>ai为突破限制)方案数+至少两个数字突破限制.........

ok,解题思路就是上面这些。但是这题比较麻烦的在于:要计算的组合数非常大且模数不为质数!!!那么Lucas定理也不管用了,我们只能用到exLucas定理。

推荐看这篇博客:https://www.cnblogs.com/candy99/p/6637629.html

直接套上大佬的exLucas会TLE,有一些地方需要优化,①Lucas函数内每次都要找MOD的因子太慢了,先预处理出来。②Fac函数内每次暴力计算阶乘太慢了,需要在调用Fac之前即C函数处先预处理阶乘。

此题代码(注释的是修改大佬模板地方,加上优化):

 #include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e5+;
ll n,m,n1,n2,MOD,a[N];
ll tot=,f[N],c[N]; ll fac[N];
void Initfac(ll p,ll pr) { //预处理阶乘
fac[]=;
for(ll i=;i<=pr;i++)
if(i%p) fac[i]=fac[i-]*i%pr;
else fac[i]=fac[i-];
} ll Pow(ll a,ll b,ll P){
ll ans=;
for(;b;b>>=,a=a*a%P)
if(b&) ans=ans*a%P;
return ans;
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y){
if(b==) d=a,x=,y=;
else exgcd(b,a%b,d,y,x),y-=(a/b)*x;
}
ll Inv(ll a,ll n){
ll d,x,y;
exgcd(a,n,d,x,y);
return d==?(x+n)%n:-;
}
ll Fac(ll n,ll p,ll pr){
if(n==) return ;
ll re=;
//for(ll i=2;i<=pr;i++) if(i%p) re=re*i%pr;
re=fac[pr];
re=Pow(re,n/pr,pr);
ll r=n%pr;
//for(int i=2;i<=r;i++) if(i%p) re=re*i%pr;
re=re*fac[r]%pr;
return re*Fac(n/p,p,pr)%pr;
}
ll C(ll n,ll m,ll p,ll pr){
if(n<m) return ;
Initfac(p,pr); //预处理阶乘
ll x=Fac(n,p,pr),y=Fac(m,p,pr),z=Fac(n-m,p,pr);
ll c=;
for(ll i=n;i;i/=p) c+=i/p;
for(ll i=m;i;i/=p) c-=i/p;
for(ll i=n-m;i;i/=p) c-=i/p;
ll a=x*Inv(y,pr)%pr*Inv(z,pr)%pr*Pow(p,c,pr)%pr;
return a*(MOD/pr)%MOD*Inv(MOD/pr,pr)%MOD;
}
ll Lucas(ll n,ll m){ //exLucas
ll x=MOD,re=;
// for(ll i=2;i<=MOD;i++) if(x%i==0){
// ll pr=1;
// while(x%i==0) x/=i,pr*=i;
// re=(re+C(n,m,i,pr))%MOD;
// }
for (int i=;i<=tot;i++) re=(re+C(n,m,f[i],c[i]))%MOD;
return re;
} void prework(ll n) {
ll t=n;
for (int i=;(ll)i*i<=n;i++) {
if (t%i==) {
tot++; f[tot]=i; c[tot]=;
while (t%i==) t/=i,c[tot]*=i;
}
}
if (t>) f[++tot]=t,c[tot]=t;
} int main()
{
int T; scanf("%d%lld",&T,&MOD);
prework(MOD);
while(T--) {
scanf("%lld%d%d%lld",&n,&n1,&n2,&m);
for(int i=;i<=n1+n2;i++) scanf("%lld",&a[i]);
for(int i=;i<=n2;i++) m-=a[n1+i]-;
m-=n; ll ans=;
for(int i=;i<(<<n1);i++) {
int tot=;
ll x=m;
for(int j=;j<n1;j++)
if (i&(<<j)) { tot++; x-=a[j+]; }
if (tot%) ans=(ans-Lucas(x+n-,n-)+MOD)%MOD;
else ans=(ans+Lucas(x+n-,n-))%MOD;
}
printf("%lld\n",ans);
}
return ;
}

记录下大佬的exLucas定理模板:

 ll Pow(ll a,ll b,ll P){
ll ans=;
for(;b;b>>=,a=a*a%P)
if(b&) ans=ans*a%P;
return ans;
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y){
if(b==) d=a,x=,y=;
else exgcd(b,a%b,d,y,x),y-=(a/b)*x;
}
ll Inv(ll a,ll n){
ll d,x,y;
exgcd(a,n,d,x,y);
return d==?(x+n)%n:-;
}
ll Fac(ll n,ll p,ll pr){
if(n==) return ;
ll re=;
for(ll i=;i<=pr;i++) if(i%p) re=re*i%pr;
re=Pow(re,n/pr,pr);
ll r=n%pr;
for(int i=;i<=r;i++) if(i%p) re=re*i%pr;
return re*Fac(n/p,p,pr)%pr;
}
ll C(ll n,ll m,ll p,ll pr){
if(n<m) return ;
ll x=Fac(n,p,pr),y=Fac(m,p,pr),z=Fac(n-m,p,pr);
ll c=;
for(ll i=n;i;i/=p) c+=i/p;
for(ll i=m;i;i/=p) c-=i/p;
for(ll i=n-m;i;i/=p) c-=i/p;
ll a=x*Inv(y,pr)%pr*Inv(z,pr)%pr*Pow(p,c,pr)%pr;
return a*(MOD/pr)%MOD*Inv(MOD/pr,pr)%MOD;
}
ll Lucas(ll n,ll m){
ll x=MOD,re=;
for(ll i=;i<=MOD;i++) if(x%i==){
ll pr=;
while(x%i==) x/=i,pr*=i;
re=(re+C(n,m,i,pr))%MOD;
}
return re;
}