J友

题目大意

你初始在(0,0)。
每个单位时间可以走向任意一个四连通格子。
问时刻n恰好在(x,y)的方案数。
答案模m。

组合数取模

显然，我们枚举向上走了多少步，就能计算出应该向左/右/下走多少步。
假设分别为p1,p2,p3,p4，那么相当于一个全排列。
由于对组合数取模比较熟悉，所以把全排列公式变成组合数相乘。
$C_{p1+p2+p3+p4}^{p1}*C_{p2+p3+p4}^{p2}*C_{p3+p4}^{p3}$
快速计算组合数取模可以参考reward题解
稍有不同的是，这题我们先枚举质因数的幂，然后预处理fac和num，再枚举向上的步数。
而且fac和num均只需要处理到第n项。
这里我们还需要多预处理一个fact[i]表示i!除去因子pp的结果。
然后做组合数的时候特判y

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll a[20],b[20],c[20],d[20],e[20],pri[32000+10],fac[1000000+10],fact[1000000+10],num[1000000+10];
bool bz[32000+10];
ll i,j,k,l,t,n,m,p,pp,x,y,top,p1,p2,p3,p4,cnt,xx,yy,ans;
ll quicksortmi(ll x,ll y,ll p){
    if (!y) return 1;
    if (y==1) return x%p;
    ll t=quicksortmi(x,y/2,p);
    t=t*t%p;
    if (y%2) t=t*(x%p)%p;
    return t;
}
void gcd(ll a,ll b){
    if (!b){
        xx=1;
        yy=0;
    }
    else{
        gcd(b,a%b);
        swap(xx,yy);
        yy-=xx*(a/b);
    }
}
ll getny(ll x,ll y){
    gcd(x,y);
    xx=(xx%y+y)%y;
    return xx;
}
ll calcfac(ll n,ll p,ll pp){
    if (n<pp) return fac[n];
    ll t=quicksortmi(fac[p-1],n/p,p);
    t=t*fac[n%p]%p;
    cnt+=n/pp;
    t=t*calcfac(n/pp,p,pp)%p;
    return t;
}
ll calc(ll x,ll y,ll p,ll pp){
    ll i,j;
    if (y<p){
        cnt=num[y]-num[x]-num[y-x];
        ll t=fact[y]*getny(fact[x]*fact[y-x]%p,p)%p*quicksortmi(pp,cnt,p)%p;
        return t;
    }
    cnt=0;
    ll A=calcfac(y,p,pp);
    ll tot=cnt;
    cnt=0;
    ll B=calcfac(x,p,pp);
    B=B*calcfac(y-x,p,pp)%p;
    B=getny(B,p);
    return A*B%p*quicksortmi(pp,tot-cnt,p)%p;
}
int main(){
    scanf("%lld%lld",&n,&m);
    scanf("%lld%lld",&x,&y);
    fo(i,2,32000){
        if (!bz[i]) pri[++k]=i;
        fo(j,1,k){
            if (pri[j]*i>32000) break;
            bz[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
    pp=m;
    fo(i,1,k){
        if (pp%pri[i]==0){
            d[++top]=1;e[top]=pri[i];
            while (pp%pri[i]==0){
                d[top]*=pri[i];
                pp/=pri[i];
            }
        }
    }
    if (pp>1){
        d[++top]=pp;
        e[top]=d[top];
    }
    fo(i,1,top) c[i]=m/d[i];
    fo(i,1,top) b[i]=getny(c[i],d[i]);
    fo(i,1,top){
        p=d[i];pp=e[i];
        fac[0]=fact[0]=1;
        num[0]=0;
        fo(j,1,min(p-1,n)){
            num[j]=0;
            k=j;
            while (k%pp==0){
                num[j]++;
                k/=pp;
            }
            num[j]+=num[j-1];
            if (j%pp==0) fac[j]=fac[j-1];
            else fac[j]=fac[j-1]*j%p;
            fact[j]=fact[j-1]*k%p;
        }
        fo(p1,0,n){
            p3=p1-x;
            if ((n-p1-p3+y)%2) continue;
            p2=(n-p1-p3+y)/2;
            p4=n-p1-p3-p2;
            if (p3<0||p2<0||p4<0) continue;
            a[i]=(a[i]+calc(p1,p1+p2+p3+p4,p,pp)*calc(p2,p2+p3+p4,p,pp)%p*calc(p3,p3+p4,p,pp)%p)%p;
        }
    }
    fo(i,1,top) ans=(ans+a[i]*b[i]%m*c[i]%m)%m;
    printf("%lld\n",ans);
}