题解（注释讲解）

#include<cstdio>
#include<algorithm>
#include<ctime>
using namespace std;const int M=190;typedef long long ll;int stat;int end;
int n;int m;ll mod;ll c[M][M];ll res;ll ans[M];ll siz;ll hsz;int a;int b;ll mi[M];
struct mar
{
 	ll mp[M][M];
 	mar(){for(int i=0;i<siz;i++)for(int j=0;j<siz;j++)mp[i][j]=0;}
 	friend mar operator *(const mar& a,const mar& b)//这里只做了4个上三角矩阵乘法，常数/2 
 	{
 		mar c;
 		for(int i=0;i<hsz;i++)//左上角 
         	for(int k=0;k<siz;k++)
 				for(int j=i;j<hsz;j++) 
 					(c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod;
 		for(int i=0;i<hsz;i++)//右上角 
         	for(int k=0;k<siz;k++)
 				for(int j=hsz+i;j<siz;j++) 
 					(c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod;
 		for(int i=hsz;i<siz;i++)//左下角 
         	for(int k=0;k<siz;k++)
 				for(int j=i-hsz;j<hsz;j++) 
 					(c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod;
 		for(int i=hsz;i<siz;i++)//左下角 
         	for(int k=0;k<siz;k++)
 				for(int j=i;j<siz;j++) 
 					(c.mp[i][j]+=a.mp[i][k]*b.mp[k][j])%=mod;
 		return c;
 	}
}tr,ret,st;
int main()
{
    scanf("%d%d%d%lld",&n,&a,&b,&mod);hsz=a+b+1;siz=2*hsz;mi[0]=1;//因为放了两个dp数组 
    for(int i=0;i<hsz;i++)c[0][i]=c[i][i]=1;//先打表组合数 
    for(int i=0;i<hsz;i++)
        for(int j=1;j<i;j++)
            c[j][i]=(c[j-1][i-1]+c[j][i-1])%mod;
    for(int i=0;i<hsz;i++)tr.mp[i][i]=1;for(int i=0;i<hsz;i++)tr.mp[hsz+i][i]=1;//单位矩阵部分 
    for(int i=0;i<hsz;i++)for(int j=i;j<hsz;j++)tr.mp[i][hsz+j]=c[i][j];//杨辉三角部分 
    for(int i=0;i<siz;i++)ret.mp[i][i]=1;
    for(int p=n;p;p>>=1,tr=tr*tr)if(p&1)ret=ret*tr;//矩阵快速幂 
    for(int i=0;i<hsz;i++)ans[i]=(ret.mp[0][i]+ret.mp[0][hsz+i])%mod;//计算y^(a+b-i)和 
    for(int i=1;i<hsz;i++){mi[i]=mi[i-1]*n%mod;}//计算幂次 
    for(int i=0;i<=a;i++)//统一使用二项式定理 
    {
        ll del=c[i][a]*mi[i]%mod*ans[a+b-i]%mod;
        if((a-i)%2){res=(res+mod-del)%mod;}else {res=(res+del)%mod;}
    }printf("%lld",res);return 0;//拜拜程序~ 
}