组合数取模(lucas定理+CRT合并)(AC)

 1 #include<bits/stdc++.h>
 2 #define re register
 3 #define int long long
 4 using namespace std;
 5 int T,n,m,mod,ans=0;
 6 inline int read(){
 7     re int a=0,b=1;re char ch=getchar();
 8     while(ch<'0'||ch>'9')
 9         b=(ch=='-')?-1:1,ch=getchar();
10     while(ch>='0'&&ch<='9')
11         a=(a<<3)+(a<<1)+(ch^48),ch=getchar();
12     return a*b;
13 }
14 inline int qpow(re int a,re int b,re int p){
15     re int res=1;a%=p;
16     for(;b;b>>=1){
17         if(b&1)res=res*a%p;
18         a=a*a%p;
19     }
20     return res%p;
21 }
22 inline int exgcd(re int a,re int b,re int &x,re int &y){
23     re int d=a;
24     if(!b){x=1,y=0;return a;}
25     d=exgcd(b,a%b,y,x);
26     y-=(a/b)*x;
27     return d;
28 }
29 inline int inv(re int a,re int m){
30     re int x,y;
31     re int d=exgcd(a,m,x,y);
32     return d==1?(m+x%m)%m:-1;
33 }
34 inline int mul(re int n,re int pi,re int pk){
35     if(!n)return 1;
36     re int ans=1;
37     if(n/pk){
38         for(re int i=2;i<=pk;i++)
39             if(i%pi) ans=ans*i%pk;
40         ans=qpow(ans,n/pk,pk);
41     }
42     for(re int i=2;i<=n%pk;i++)
43         if(i%pi) ans=ans*i%pk;
44     return ans*mul(n/pi,pi,pk)%pk;
45 }
46 inline int C(re int n,re int m,re int p,re int pi,re int pk)
47 {
48     if(m>n)return 0;
49     re int a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk);
50     re int k=0,ans;
51     for(re int i=n;i;i/=pi) k+=i/pi;
52     for(re int i=m;i;i/=pi) k-=i/pi;
53     for(re int i=n-m;i;i/=pi) k-=i/pi;
54     ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qpow(pi,k,pk)%pk;
55     ans=ans*(p/pk)%p*inv(p/pk,pk)%p;
56     return ans;
57 }
58 inline int lucas(re int m,re int n,re int p){
59     re int x=p,ans=0;
60     for(re int i=2;i<=p;i++){
61         if(!(x%i)){
62             re int pk=1;
63             while(!(x%i))pk*=i,x/=i;
64             ans=(ans+C(n,m,p,i,pk))%p;
65         }
66     }
67     return ans;
68 }
69 signed main(){
70     T=read(),mod=read();
71     n=read(),m=read();
72     if(mod==1||T<n+m){puts("0");return 0;}
73     for(re int k=m;k<=(T+m-n)/2;k++)
74         ans+=lucas(k,T,mod)%mod*lucas(k-m,T-k,mod)%mod*lucas((T+m+n-(k<<1))/2,T-(k<<1)+m,mod)%mod;
75     printf("%lld\n",ans%mod);
76     return 0;
77 }
 1 #include<bits/stdc++.h>
 2 #define re register
 3 #define ll long long
 4 using namespace std;
 5 ll T,n,m,mod,ans=0;
 6 inline ll read(){
 7     re ll a=0,b=1;re char ch=getchar();
 8     while(ch<'0'||ch>'9')
 9         b=(ch=='-')?-1:1,ch=getchar();
10     while(ch>='0'&&ch<='9')
11         a=(a<<3)+(a<<1)+(ch^48),ch=getchar();
12     return a*b;
13 }
14 inline void exgcd(re ll &x,re ll &y,re ll a,re ll b)
15 {
16     if(!b){x=1,y=0;return;}
17     exgcd(x,y,b,a%b);
18     ll tmp=x;x=y;
19     y=tmp-a/b*y;
20 }
21 inline ll inv(re ll a,re ll p)
22 {
23     re ll x,y;
24     exgcd(x,y,a,p);
25     return x;
26 }
27 inline ll qpow(re ll a,re ll b,re ll p){
28     re ll res=1;a%=p;
29     for(;b;b>>=1){
30         if(b&1)res=res*a%p;
31         a=a*a%p;
32     }
33     return res%p;
34 }
35 inline ll fac(re ll n,re ll p,re ll tp)
36 {
37     if(!n) return 1;
38     re ll f=1,res=1;
39     for(re ll i=1;i<tp;i++)
40     {
41         if(i%p) (f*=i)%=tp;
42         if(i==n%tp) res=f;
43     }
44     f=qpow(f,n/tp,tp);
45     return fac(n/p,p,tp)*f%tp*res%tp;
46 }
47 inline ll C(re ll m,re ll n,re ll p,re ll tp)
48 {
49     re ll ct=0;
50     for(re ll i=m;i;i/=p) ct+=i/p;
51     for(re ll i=n;i;i/=p) ct-=i/p;
52     for(re ll i=m-n;i;i/=p) ct-=i/p;
53     return fac(m,p,tp)*inv(fac(n,p,tp),tp)%tp*inv(fac(m-n,p,tp),tp)%tp*qpow(p,ct,tp)%tp;
54 }
55 inline ll CRT(re ll m,re ll n,re ll p,re ll tp)
56 {
57     return C(m,n,p,tp)*(mod/tp)%mod*inv(mod/tp,tp)%mod;
58 }
59 inline ll lucas(re ll n,re ll m,re ll p)
60 {
61     ll ans=0;
62     for(re ll i=2;i*i<=p;i++)
63     {
64         if(p%i==0)
65         {
66             re ll tp=1;
67             while(p%i==0) tp*=i,p/=i;
68             (ans+=CRT(m,n,i,tp))%=mod;
69         }
70     }
71     if(p!=1) (ans+=CRT(m,n,p,p))%=mod;
72     return (ans%mod+mod)%mod;
73 }
74 signed main(){
75     T=read(),mod=read();
76     n=read(),m=read();
77     if(mod==1||T<n+m){puts("0");return 0;}
78     for(re int k=m;k<=(T+m-n)/2;k++)
79         ans+=lucas(k,T,mod)%mod*lucas(k-m,T-k,mod)%mod*lucas((T+m+n-(k<<1))/2,T-(k<<1)+m,mod)%mod;
80     printf("%lld\n",ans%mod);
81     return 0;
82 }

 以下是ac代码

附上新的公式推导方式(经实测可AC)

设k为向上走的步数(纯的,不包括向下),那么向下走的步数为k-m

即C(T,k)×C(T-k,k-m);

由上式可推出 向右的与向左的步数之和为T-2*k+m,而向右的比向左的多n,利用和差公式,易得,向右的为(T-2*k+m+n)/2;

且k>=m,(T-2*k+m+n)/2>=n;可得,k->[m,(T+m-n)/2];

故整体公式为

  ∑(k->[m,(T+m-n)/2])(C(T,k)×C(T-k,k-m)×C(T-2*k+m,(T-2*k+m+n)/2));

 1 #include <cmath>
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <iostream>
 5 using namespace std;
 6 #define QAQ 100010
 7 #define re register
 8 #define int long long
 9 #define ll long long
10 #define fup(i,a,b) for(re int i=a;i<=b;++i)
11 #define fdn(i,a,b) for(re int i=a;i>=b;--i)
12 int n,m,mod,t;
13 ll ans;
14 int fac[QAQ],inv[QAQ];
15 int tot,p[QAQ],w[QAQ];
16 inline int quick_fk_me(re int a,re int b,re int pk){
17     int ans=1;
18     for(;b;b>>=1){
19         if(b&1)ans=1ll*ans*a%pk;
20         a=1ll*a*a%pk;
21     }
22 }
23 void exgcd(re ll a,re ll b,re ll &x,re ll &y){
24     if(b==0){x=1,y=0;return;}
25     exgcd(b,a%b,x,y);
26     re int z=x;
27     x=y,y=z-a/b*y;
28 }
29 inline ll crt(){
30     re ll ans=0,lcm=1,x,y;
31     fup(i,1,tot)lcm*=p[i];
32     fup(i,1,tot){
33         re ll tmp=lcm/p[i];
34         exgcd(tmp,p[i],x,y);
35         x=(x%p[i]+p[i])%p[i];
36         ans=(ans+tmp*x*w[i])%lcm;
37     }
38     return (ans+lcm)%lcm;
39 }
40 inline void Get_fac_and_inv(re int top,re int pk){
41     fac[0]=fac[1]=inv[0]=inv[1]=1;
42     fup(i,2,top)fac[i]=1ll*fac[i-1]*i%pk;
43     fup(i,2,top)inv[i]=(pk-pk/i)*inv[pk%i]%pk;
44     fup(i,1,top)inv[i]=1ll*inv[i-1]*inv[i]%pk;
45 }
46 int CC(re int n,re int m,re int pk){
47     if(m>n)return 0;
48     return fac[n]*inv[m]%pk*inv[n-m]%pk;
49 }
50 int lucas(re int n,re int m,re int pk){
51     if(!m)return 1;
52     return CC(n%pk,m%pk,pk)*lucas(n/pk,m/pk,pk)%pk;
53 }
54 inline bool divide(int n){
55     re int top=sqrt(n);
56     fup(i,2,top){
57         if(n%i)continue;
58         p[++tot]=i;
59         while(n%i==0)n/=i;
60     }
61     if(n>1)p[++tot]=n;
62     return tot==1;
63 }
64 main(){
65     scanf("%lld%lld",&t,&mod);
66     scanf("%lld%lld",&n,&m);
67     if(n<0)n=-n;
68     if(m<0)m=-m;
69     re int st=m,ed=(t+m-n)>>1;
70     if(divide(mod)){
71         Get_fac_and_inv(min(t,mod),mod);//什么纱吊错 妈耶 不求阶乘??!!
72         fup(k,st,ed)
73         ans=(ans+lucas(t,k,mod)*lucas(t-k,k-m,mod)%mod*lucas(t-2*k+m,(t-2*k+m-n)>>1,mod)%mod)%mod;
74         printf("%lld\n",ans);
75         return 0;
76     }
77     fup(i,1,tot){
78         re int pk=p[i];
79         Get_fac_and_inv(min(t,pk),pk);
80         fup(k,st,ed)
81         w[i]=(w[i]+lucas(t,k,pk)*lucas(t-k,k-m,pk)%pk*lucas(t-2*k+m,(t-2*k+m-n)>>1,pk)%pk)%pk;
82     }
83     printf("%lld\n",crt());
84 }
 1 #include<bits/stdc++.h>
 2 #define re register
 3 #define int long long
 4 using namespace std;
 5 int T,n,m,mod,ans=0,tot=0,fac[100001],pi,pk,p;
 6 inline int read(){
 7     re int a=0,b=1;re char ch=getchar();
 8     while(ch<'0'||ch>'9')
 9         b=(ch=='-')?-1:1,ch=getchar();
10     while(ch>='0'&&ch<='9')
11         a=(a<<3)+(a<<1)+(ch^48),ch=getchar();
12     return a*b;
13 }
14 inline int qpow(re int a,re int b){
15     re int res=1;a%=pk;
16     for(;b;b>>=1){
17         if(b&1)res=res*a%pk;
18         a=a*a%pk;
19     }
20     return res%pk;
21 }
22 inline int getcnt(re int x){
23     return x?x/pi+getcnt(x/pi):0;
24 }
25 inline int getsum(re int x){
26     return x?((x/pk?qpow(fac[pk],x/pk):1)*fac[x%pk]%pk*getsum(x/pi)%pk):1;
27 } 
28 inline int C(re int x,re int y){
29     if(y>x)return 0;
30     re int a0=getsum(x),b0=getcnt(x);
31     re int a1=getsum(y),b1=getcnt(y);
32     re int a2=getsum(x-y),b2=getcnt(x-y);
33       b0-=b1+b2;
34       if(b0>=tot)return 0;//含质因子的个数
35       return a0*qpow(a1,p-1)%pk*qpow(a2,p-1)%pk*qpow(pi,b0)%pk;
36 }
37 inline int work(){
38     re int res=0,a,b,c,d;
39       fac[0]=1;
40       for(re int i=1;i<=pk&&i<=T;i++)
41         fac[i]=(i%pi)?fac[i-1]*i%pk:fac[i-1];//处理阶乘
42       a=0,b=n,c=(T-n-m)>>1,d=(T+m-n)>>1;//各方向初值
43       while(c>=0){
44         res=(res+C(T,a+b)*C(a+b,a)%pk*C(c+d,c)%pk)%pk;
45         ++a,++b,--c,--d;
46      }//处理组合数,上下左右的道路可能性
47       return res;
48 }
49 signed main(){
50     T=read(),mod=read(),n=read(),m=read();
51     n=(n<0)?-n:n;m=(m<0)?-m:m;//处理负数
52     re int l=mod;
53     for(re int i=2;l>1;i++){
54         if(i*i>l)i=l;//分解质因数
55         if(l%i==0){
56             pi=i,pk=1,tot=0;
57             while(l%pi==0)
58                 l/=pi,pk*=pi,++tot;
59             p=pk/pi*(pi-1);
60             (ans+=work()%mod*(mod/pk)%mod*qpow(mod/pk,p-1)%mod)%=mod;//CRT
61         }
62     }
63     printf("%lld\n",ans);
64     return 0;
65 }

另有 学习笔记“组合数取模”推荐

转载于:https://www.cnblogs.com/Hzoi-lyl/p/11226099.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值