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

题意:给定方程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 }
View Code

记录下大佬的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 }
View Code
原文地址:https://www.cnblogs.com/clno1/p/11662402.html