题目:http://www.51nod.com/Challenge/Problem.html#!#problemId=1362
首先,\( f[i][j] \) 是一个 \( i \) 次多项式;
如果考虑差分,用一个列向量维护 0 次差分到 \( n \) 次差分即可,在第 \( n \) 次上差分数组已经是一个常数;
最后一行再维护一个 0 次差分的前缀和,0 次差分其实就是答案;
为了预处理 0 位置上的各次差分值,一开始先 n^2 求出 \( f[0][0] \) 到 \( f[n][n] \),差分一下得到初始矩阵;
转移就是本层加上下一层的差分值,得到本层的下一个位置,\( n \) 次的差分值不变;
需要注意快速幂时,子函数是不能传超过大约 500*500 的数组的,所以把矩阵开成全局的;
可惜这样是 \( n^{3}logm \),会TLE;
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=805; int n,m,f[xn][xn],d[xn][xn],P; int upt(int x){while(x>=P)x-=P; while(x<0)x+=P; return x;} struct N{int a[xn][xn];N(){memset(a,0,sizeof a);}void init(){for(int i=0;i<=n+1;i++)a[i][i]=1;}void clr(){memset(a,0,sizeof a);}N operator * (const N &y) const{N ret;for(int i=0;i<=n+1;i++)for(int k=0;k<=n+1;k++)for(int j=0;j<=n+1;j++)ret.a[i][j]=upt(ret.a[i][j]+(ll)a[i][k]*y.a[k][j]%P);return ret;} }a,g,t; void print(N a) {for(int i=0;i<=n+1;i++,puts(""))for(int j=0;j<=n+1;j++)printf("%d",a.a[i][j]); } void pw(int b)// {t.init();for(;b;b>>=1,g=g*g)if(b&1)t=t*g; } int main() {while(scanf("%d%d%d",&n,&m,&P)!=EOF){memset(f,0,sizeof f);memset(d,0,sizeof d);f[0][0]=1;for(int i=0;i<=n;i++)for(int j=0;j<=n;j++){if(i)f[i][j]=upt(f[i][j]+f[i-1][j]);if(j)f[i][j]=upt(f[i][j]+f[i][j-1]);if(i&&j)f[i][j]=upt(f[i][j]+f[i-1][j-1]);}for(int j=0;j<=n;j++)d[0][j]=f[n][j];for(int i=1;i<=n;i++)for(int j=0;j<=n-i;j++)d[i][j]=upt(d[i-1][j+1]-d[i-1][j]);a.clr(); g.clr(); t.clr();for(int i=0;i<=n;i++)a.a[i][0]=d[i][0];for(int i=0;i<n;i++)g.a[i][i]=g.a[i][i+1]=1;g.a[n][n]=1;g.a[n+1][0]=g.a[n+1][n+1]=1;pw(m);int sum=0;for(int i=0;i<=n+1;i++)sum=upt(sum+(ll)t.a[0][i]*a.a[i][0]%P),sum=upt(sum+(ll)t.a[n+1][i]*a.a[i][0]%P);printf("%d\n",sum);}return 0; }
其实可以考虑组合数:
因为斜着走既向下又向右,不好判断,所以不妨枚举斜着走了几格,假设 \( n<=m \),得到
\( ans = \sum\limits_{j=0}^{m} \sum\limits_{i=0}^{n} C_{i+j}^{i} * C_{j}^{n-i} \)
即 \( ans = \sum\limits_{j=0}^{m} \sum\limits_{i=0}^{n} C_{n}^{i} * C_{i+j}^{n} \),或者其实可以直接写出这个式子;
然后 \( ans = \sum\limits_{i=0}^{n} C_{n}^{i} * \sum\limits_{j=0}^{m} C_{i+j}^{n} \)
\( ans = \sum\limits_{i=0}^{n} C_{n}^{i} * C_{i+m+1}^{n+1} \)
于是求组合数即可;
但模数不是质数,没有逆元;
可以像扩展 Lucas 的做法一样,提取出模数的质因子,剩余的部分就和模数互质,可以用 exgcd 求逆元;
质因子的部分直接把次数加减即可。
代码如下:
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=805,xm=35; int n,m,mod,p[xm],cnt,t[xm]; void div(int x) {for(int i=2;i*i<=x;i++){if(x%i)continue;p[++cnt]=i;while(x%i==0)x/=i;}if(x>1)p[++cnt]=x; } ll pw(ll a,int b) {ll ret=1;for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod;return ret; } void exgcd(int a,int b,int &x,int &y) {if(!b){x=1; y=0; return;}exgcd(b,a%b,y,x);y-=a/b*x; } int inv(int a,int b){int x,y; exgcd(a,b,x,y); return (x%b+b)%b;}//%b int get(int x,int tp) {for(int i=1;i<=cnt;i++){if(x%p[i])continue;int cnt=0;while(x%p[i]==0)cnt++,x/=p[i];t[i]+=cnt*tp;}return x; } int C(int n,int m) {if(n<m)return 0;//!!if(m==0)return 1;memset(t,0,sizeof t);int ret=1;for(int i=n-m+1;i<=n;i++)ret=(ll)ret*get(i,1)%mod;for(int i=1;i<=m;i++)ret=(ll)ret*inv(get(i,-1),mod)%mod;for(int i=1;i<=cnt;i++)ret=(ll)ret*pw(p[i],t[i])%mod;return ret; } int main() {while(scanf("%d%d%d",&n,&m,&mod)==3){cnt=0; div(mod); int ans=0;for(int i=0;i<=n;i++)ans=(ans+(ll)C(n,i)*C(i+m+1,n+1))%mod;printf("%d\n",ans);}return 0; }