题目：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 ",sum); } return 0; }
其实可以考虑组合数：
因为斜着走既向下又向右，不好判断，所以不妨枚举斜着走了几格，假设 ( n<=m )，得到
( ans = sumlimits_{j=0}^{m} sumlimits_{i=0}^{n} C_{i+j}^{i} * C_{j}^{n-i} )
即 ( ans = sumlimits_{j=0}^{m} sumlimits_{i=0}^{n} C_{n}^{i} * C_{i+j}^{n} )，或者其实可以直接写出这个式子；
然后 ( ans = sumlimits_{i=0}^{n} C_{n}^{i} * sumlimits_{j=0}^{m} C_{i+j}^{n} )
( ans = sumlimits_{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 ",ans); } return 0; }