HZOJ visit

最新推荐文章于 2021-05-24 16:36:19 发布

weixin_30415801

最新推荐文章于 2021-05-24 16:36:19 发布

阅读量93

点赞数

原文链接：http://www.cnblogs.com/Al-Ca/p/11228419.html

版权

本文介绍了一种解决复杂路径计数问题的算法，分为三个部分：对于前30%的数据采用动态规划方法；对于中间30%的数据通过组合数学简化计算；对于剩余的数据集，当模数不是质数时，利用中国剩余定理解决组合数计算问题。

摘要由CSDN通过智能技术生成

对于前30%的数据，可以考虑dp，f[i][j][k]表示时间为i，在i，j位置的方案数，枚举转移即可。要注意的是可以走到矩阵外。

对于另外30%数据，考虑推一下式子，设向右走y步，左z，上s，下x。那么y-z=n,s-x=m。所以我们枚举s就可以求得sxzy，步数确定之后就比较简单了，显然答案为

$∑C_T^s*C_{T-s}^x*C_{T-s-x}^z*C_{T-s-x-z}^y$,化减得ans=∑$\frac{T!}{s!x!z!y!}$,由于mod是质数，直接逆元干就行了。

对于100%数据，难点就在于p不是质数，开始我想着分解质因数，然后T了。

考虑一下如何求组合数%合数？（合数为若干质数乘积）

将合数质因数分解，用卢卡斯求出组合数mod每个质数，然后发现其实就是线性同余方程组，CRT解即可。

  1 #include<iostream>
  2 #include<cstring>
  3 #include<cstdio>
  4 #include<queue>
  5 #define int LL
  6 #define ma(x) memset(x,0,sizeof(x))
  7 #define LL long long
  8 using namespace std;
  9 int T,mod,n,m;
 10 int s,x,z,y;
 11 LL f[210][310][310];
 12 LL jc[100010];
 13 int cnt[100010];
 14 int prime[100010],num;
 15 bool isprime[100010];
 16 LL exgcd(int a,int b,int &x,int &y)
 17 {
 18     if(!b){x=1,y=0;return a;}    
 19     int gcd=exgcd(b,a%b,x,y),t=x;
 20     x=y,y=t-a/b*y;
 21     return gcd;
 22 }
 23 void ss()
 24 {
 25     const int N=100010;
 26     for(int i=2;i<=N;i++)isprime[i]=1;
 27     for(int i=2;i<=N;i++)    
 28     {
 29         if(isprime[i])prime[++num]=i;
 30         for(int j=1;j<=num&&i*prime[j]<=N;j++)
 31         {
 32             isprime[i*prime[j]]=0;
 33             if(i%prime[j]==0)break;    
 34         }
 35     }
 36 }
 37 void add(int x,int nu)
 38 {
 39     for(int i=1;prime[i]*prime[i]<=x;i++)
 40     while(x%prime[i]==0){cnt[prime[i]]+=nu;x/=prime[i];}
 41     cnt[x]+=nu;
 42 }
 43 bool pdprime(int x)
 44 {
 45     for(int i=2;i*i<=x;i++)
 46         if(x%i==0)return 0;
 47     return 1;
 48 }
 49 int fj[1000010],cntt;
 50 LL poww(LL a,int b,int mod)
 51 {
 52     LL ans=1;
 53     while(b)
 54     {
 55         if(b&1)ans=ans*a%mod;
 56         a=a*a%mod;
 57         b=b>>1;
 58     }
 59     return ans;
 60 }
 61 LL C(int n,int m,int mod)
 62 {
 63     if(m>n)return 0;
 64     return jc[n]*poww(jc[m],mod-2,mod)%mod*poww(jc[n-m],mod-2,mod)%mod;
 65 }
 66 LL Lucas(int n,int m,int mod)
 67 {
 68     if(m>n)return 0;
 69     if(!m)return 1;
 70     return Lucas(n/mod,m/mod,mod)*C(n%mod,m%mod,mod)%mod;
 71 }
 72 int CRT(int W[],int B[],int k)
 73 {
 74     int x,y,a=0,m,n=1;
 75     for(int i=1;i<=k;i++)
 76         n*=W[i];
 77     for(int i=1;i<=k;i++)
 78     {
 79         m=n/W[i];
 80         exgcd(W[i],m,x,y);
 81         a=(a+y*m*B[i])%n;
 82     }
 83     return a>0?a:(a+n);
 84 }
 85 int w[10000],b[10000];
 86 signed main()
 87 {
 88 //    freopen("out.doc","w",stdout);
 89 
 90     ss();
 91     cin>>T>>mod>>n>>m;int tem=mod;
 92     for(int i=1;prime[i]*prime[i]<=tem;i++)
 93     while(tem%prime[i]==0)
 94     {
 95         fj[++cntt]=prime[i];
 96         tem/=prime[i];
 97     }
 98     if(tem>1)fj[++cntt]=tem;
 99 //    for(int i=1;i<=cntt;i++)cout<<fj[i]<<" ";puts("");
100     if(T<=100)
101     {
102         f[0][100][100]=1;
103         for(int i=1;i<=T;i++)
104             for(int j=1;j<=n+200;j++)
105                 for(int k=1;k<=m+200;k++)
106                     f[i][j][k]=((f[i-1][j][k-1]+f[i-1][j][k+1])%mod+(f[i-1][j-1][k]+f[i-1][j+1][k])%mod)%mod;
107         cout<<f[T][n+100][m+100]%mod<<endl;
108         return 0;
109     }
110     if(pdprime(mod))
111     {
112         if(n<0)n=-n;
113         if(m<0)m=-m;
114         jc[0]=1;for(int i=1;i<=100000;i++)jc[i]=jc[i-1]*i%mod;
115         LL ans=0;
116         for(int s=m;s<=T;s++)
117         {
118             x=s-m;
119             if((T-s-x+n)%2!=0)continue;
120             y=(T-s-x+n)/2;
121             z=y-n;
122             if(y<0||z<0)break;
123             ans=(ans+jc[T]*poww(jc[s],mod-2,mod)%mod*poww(jc[x],mod-2,mod)%mod*poww(jc[z],mod-2,mod)%mod*poww(jc[y],mod-2,mod)%mod)%mod;
124         }
125         cout<<ans<<endl;
126     }
127     else
128     {
129         if(n<0)n=-n;
130         if(m<0)m=-m;
131         jc[0]=1;for(int i=1;i<=100000;i++)jc[i]=jc[i-1]*i%mod;
132         LL ans=0;
133         for(int s=m;s<=T;s++)
134         {
135             x=s-m;
136             if((T-s-x+n)%2!=0)continue;
137             y=(T-s-x+n)/2;
138             z=y-n;
139             if(y<0||z<0)break;
140             LL tem=1;
141             ma(w),ma(b);
142             for(int i=1;i<=cntt;i++)    
143                 w[i]=fj[i],b[i]=Lucas(T,s,fj[i]);
144             tem=tem*CRT(w,b,cntt)%mod;
145             ma(w),ma(b);
146             for(int i=1;i<=cntt;i++)    
147                 w[i]=fj[i],b[i]=Lucas(T-s,x,fj[i]);
148             tem=tem*CRT(w,b,cntt)%mod;
149             ma(w),ma(b);
150             for(int i=1;i<=cntt;i++)    
151                 w[i]=fj[i],b[i]=Lucas(T-s-x,z,fj[i]);
152             tem=tem*CRT(w,b,cntt)%mod;
153             ans=(ans+tem)%mod;
154         }
155         cout<<ans%mod<<endl;
156     }
157 }