题意:给定方程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函数处先预处理阶乘。
此题代码(注释的是修改大佬模板地方,加上优化):
1 #include <bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const int N=1e5+10; 5 ll n,m,n1,n2,MOD,a[N]; 6 ll tot=0,f[N],c[N]; 7 8 ll fac[N]; 9 void Initfac(ll p,ll pr) { //预处理阶乘 10 fac[0]=1; 11 for(ll i=1;i<=pr;i++) 12 if(i%p) fac[i]=fac[i-1]*i%pr; 13 else fac[i]=fac[i-1]; 14 } 15 16 ll Pow(ll a,ll b,ll P){ 17 ll ans=1; 18 for(;b;b>>=1,a=a*a%P) 19 if(b&1) ans=ans*a%P; 20 return ans; 21 } 22 void exgcd(ll a,ll b,ll &d,ll &x,ll &y){ 23 if(b==0) d=a,x=1,y=0; 24 else exgcd(b,a%b,d,y,x),y-=(a/b)*x; 25 } 26 ll Inv(ll a,ll n){ 27 ll d,x,y; 28 exgcd(a,n,d,x,y); 29 return d==1?(x+n)%n:-1; 30 } 31 ll Fac(ll n,ll p,ll pr){ 32 if(n==0) return 1; 33 ll re=1; 34 //for(ll i=2;i<=pr;i++) if(i%p) re=re*i%pr; 35 re=fac[pr]; 36 re=Pow(re,n/pr,pr); 37 ll r=n%pr; 38 //for(int i=2;i<=r;i++) if(i%p) re=re*i%pr; 39 re=re*fac[r]%pr; 40 return re*Fac(n/p,p,pr)%pr; 41 } 42 ll C(ll n,ll m,ll p,ll pr){ 43 if(n<m) return 0; 44 Initfac(p,pr); //预处理阶乘 45 ll x=Fac(n,p,pr),y=Fac(m,p,pr),z=Fac(n-m,p,pr); 46 ll c=0; 47 for(ll i=n;i;i/=p) c+=i/p; 48 for(ll i=m;i;i/=p) c-=i/p; 49 for(ll i=n-m;i;i/=p) c-=i/p; 50 ll a=x*Inv(y,pr)%pr*Inv(z,pr)%pr*Pow(p,c,pr)%pr; 51 return a*(MOD/pr)%MOD*Inv(MOD/pr,pr)%MOD; 52 } 53 ll Lucas(ll n,ll m){ //exLucas 54 ll x=MOD,re=0; 55 // for(ll i=2;i<=MOD;i++) if(x%i==0){ 56 // ll pr=1; 57 // while(x%i==0) x/=i,pr*=i; 58 // re=(re+C(n,m,i,pr))%MOD; 59 // } 60 for (int i=1;i<=tot;i++) re=(re+C(n,m,f[i],c[i]))%MOD; 61 return re; 62 } 63 64 void prework(ll n) { 65 ll t=n; 66 for (int i=2;(ll)i*i<=n;i++) { 67 if (t%i==0) { 68 tot++; f[tot]=i; c[tot]=1; 69 while (t%i==0) t/=i,c[tot]*=i; 70 } 71 } 72 if (t>1) f[++tot]=t,c[tot]=t; 73 } 74 75 int main() 76 { 77 int T; scanf("%d%lld",&T,&MOD); 78 prework(MOD); 79 while(T--) { 80 scanf("%lld%d%d%lld",&n,&n1,&n2,&m); 81 for(int i=1;i<=n1+n2;i++) scanf("%lld",&a[i]); 82 for(int i=1;i<=n2;i++) m-=a[n1+i]-1; 83 m-=n; 84 85 ll ans=0; 86 for(int i=0;i<(1<<n1);i++) { 87 int tot=0; 88 ll x=m; 89 for(int j=0;j<n1;j++) 90 if (i&(1<<j)) { tot++; x-=a[j+1]; } 91 if (tot%2) ans=(ans-Lucas(x+n-1,n-1)+MOD)%MOD; 92 else ans=(ans+Lucas(x+n-1,n-1))%MOD; 93 } 94 printf("%lld ",ans); 95 } 96 return 0; 97 }
记录下大佬的exLucas定理模板:
1 ll Pow(ll a,ll b,ll P){ 2 ll ans=1; 3 for(;b;b>>=1,a=a*a%P) 4 if(b&1) ans=ans*a%P; 5 return ans; 6 } 7 void exgcd(ll a,ll b,ll &d,ll &x,ll &y){ 8 if(b==0) d=a,x=1,y=0; 9 else exgcd(b,a%b,d,y,x),y-=(a/b)*x; 10 } 11 ll Inv(ll a,ll n){ 12 ll d,x,y; 13 exgcd(a,n,d,x,y); 14 return d==1?(x+n)%n:-1; 15 } 16 ll Fac(ll n,ll p,ll pr){ 17 if(n==0) return 1; 18 ll re=1; 19 for(ll i=2;i<=pr;i++) if(i%p) re=re*i%pr; 20 re=Pow(re,n/pr,pr); 21 ll r=n%pr; 22 for(int i=2;i<=r;i++) if(i%p) re=re*i%pr; 23 return re*Fac(n/p,p,pr)%pr; 24 } 25 ll C(ll n,ll m,ll p,ll pr){ 26 if(n<m) return 0; 27 ll x=Fac(n,p,pr),y=Fac(m,p,pr),z=Fac(n-m,p,pr); 28 ll c=0; 29 for(ll i=n;i;i/=p) c+=i/p; 30 for(ll i=m;i;i/=p) c-=i/p; 31 for(ll i=n-m;i;i/=p) c-=i/p; 32 ll a=x*Inv(y,pr)%pr*Inv(z,pr)%pr*Pow(p,c,pr)%pr; 33 return a*(MOD/pr)%MOD*Inv(MOD/pr,pr)%MOD; 34 } 35 ll Lucas(ll n,ll m){ 36 ll x=MOD,re=0; 37 for(ll i=2;i<=MOD;i++) if(x%i==0){ 38 ll pr=1; 39 while(x%i==0) x/=i,pr*=i; 40 re=(re+C(n,m,i,pr))%MOD; 41 } 42 return re; 43 }