题目：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; }
