HZOJ visit

浅浅的花香味﹌ 2021-10-26 11:44 294阅读 0赞

对于前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了。

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

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

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

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 }

转载于:https://www.cnblogs.com/Al-Ca/p/11228419.html

[Link 1]: https://blog.csdn.net/qq_40772692/article/details/81902590
[ContractedBlock.gif]: https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif
[ExpandedBlockStart.gif]: /images/20211026/945f12707783444c90784b32ebd092c1.png