[51nod1362]搬箱子

Description

一个n*m的棋盘，每一步可以从(x,y)走到(x,y+1)或(x+1,y)或(x+1,y+1).
求从(0,0)走到最后一行的方案数，答案对p取模。
n<=800，m,p<=10^9

Solution

显然可以枚举斜走的步数。
然后再枚举走到(n,j)，我们要有梦想这个一定是可以化简的。
那么

$$Ans=\sum_{i=0}^{min(n,m)}\sum_{j=i}^{m}C_{n-2i+j}^{n-i}*C_{n-i+j}^{i}$$
把后面的东西回归本源 $$Ans=\sum_{i=0}^{min(n,m)}\sum_{j=i}^{m}{{(n-2i+j)!(n-i+j)!}\over {(n-i)!(j-i)!i!(n-2i+j)!}}$$
发现可以约掉，然后套用 $${^(a+b+c)!\over {a!b!c!}}=C_{a+b}^{a}*C_{a+b+c}^{a+b}$$
就可以得到 $$Ans=\sum_{i=0}^{min(n,m)}C_n^{i}\sum_{j=i}^{m}C_{n-i+j}^{n}$$
然后因为 $\sum_{i=0}^{n}C_{i}^{k}=C_{n+1}^{k+1}$
所以 $$Ans=\sum_{i=0}^{min(n,m)}C_{n}^{i}C_{n+m-i+1}^{n+1}$$
现在的问题就是如何计算这个组合数。
模数不一定是质数！
我们发现前一个可以直接预处理，后一个只有n+1项。
那么我们可以维护所有阶乘中出现过的p的质因数的次数，然后用上面的减去下面的。
剩下的直接算就好了。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=1605;
int n,m,p,ans,a[N],b[N],c[N][N],d[N][N],f[N],g[N];
int mi(int x,int y) {
    int z=1;
    for(;y;y/=2,x=(ll)x*x%p) if (y&1) z=(ll)z*x%p;
    return z;
}
void exgcd(int a,int b,int &x,int &y) {
    if (!b) {
        x=1;y=0;
        return;
    }
    int xx,yy;
    exgcd(b,a%b,xx,yy);
    x=yy;y=xx-a/b*yy;
}
int inv(int x) {
    int a,b;
    exgcd(x,p,a,b);
    ((a%=p)+=p)%=p;
    return a;
}
int main() {
    scanf("%d%d%d",&n,&m,&p);
    if (n>=m) {
        c[0][0]=1;
        fo(i,1,n+m+1) {
            c[i][0]=1;
            fo(j,1,i) c[i][j]=(c[i-1][j]+c[i-1][j-1])%p;
        }
        fo(i,0,min(n,m)) (ans+=(ll)c[n][i]*c[n+m-i+1][n+1]%p)%=p;
        printf("%d\n",ans);
        return 0;
    }
    c[0][0]=1;
    fo(i,1,n) {
        c[i][0]=1;
        fo(j,1,i) c[i][j]=(c[i-1][j]+c[i-1][j-1])%p;
    }
    int k=p;f[0]=g[0]=1;
    fo(i,2,sqrt(p)) 
        if (!(k%i)) {
            a[++a[0]]=i;
            while (!(k%i)) k/=i;
        }
    if (k>1) a[++a[0]]=k;
    fo(i,1,n+1) {
        int x=i;
        fo(j,1,a[0]) 
            while (!(x%a[j])) x/=a[j],b[j]++;
        f[i]=(ll)f[i-1]*inv(x)%p;
    }
    int l=m-n+1,r=m+n+1;
    fo(i,l,r) {
        int x=i;
        fo(j,1,a[0]) 
            while (!(x%a[j])) x/=a[j],d[i-l+1][j]++;
        g[i-l+1]=(ll)g[i-l]*x%p;
    }
    fo(i,1,r-l+1) fo(j,1,a[0]) d[i][j]+=d[i-1][j];
    fo(i,0,n) {
        int sum=c[n][i],res=1;
        sum=(ll)sum*f[n+1]%p*g[m+n+1-i-l+1]%p*inv(g[m-i-l+1])%p;
        fo(j,1,a[0]) {
            int k=-b[j];
            k+=d[m+n+1-i-l+1][j]-d[m-i-l+1][j];
            res=(ll)res*mi(a[j],k)%p;
        }
        (ans+=(ll)sum*res%p)%=p;
    }
    printf("%d\n",ans);
}