解题思路:
xayb=(n−y)ayb=∑i=0a(ai)ni(−1)a−iya+b−i
x
a
y
b
=
(
n
−
y
)
a
y
b
=
∑
i
=
0
a
(
a
i
)
n
i
(
−
1
)
a
−
i
y
a
+
b
−
i
设
f[k][i][0/1]
f
[
k
]
[
i
]
[
0
/
1
]
表示长度为
k
k
,结尾为的序列中所有合法的
y
y
的次方之和,则:
f[k][i][0]=f[k−1][i][0]+f[k−1][i][1]
f
[
k
]
[
i
]
[
0
]
=
f
[
k
−
1
]
[
i
]
[
0
]
+
f
[
k
−
1
]
[
i
]
[
1
]
f[k][i][1]=∑(yl+1)i=∑j=0if[k−1][j][0]
f
[
k
]
[
i
]
[
1
]
=
∑
(
y
l
+
1
)
i
=
∑
j
=
0
i
f
[
k
−
1
]
[
j
]
[
0
]
直接矩阵快速幂算出后统计答案即可。
#include<bits/stdc++.h>
#define ll unsigned long long
using namespace std;
int getint()
{
int i=0,f=1;char c;
for(c=getchar();(c!='-')&&(c<'0'||c>'9');c=getchar());
if(c=='-')c=getchar(),f=-1;
for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
return i*f;
}
const int N=200;
ll n,m,a,b,mod,ans,c[N][N];
inline void add(ll &x,ll y){x=x+y>=mod?x+y-mod:x+y;}
struct matrix
{
ll a[N][N];
matrix(){memset(a,0,sizeof(a));}
inline friend matrix operator * (const matrix &A,const matrix &B)
{
matrix res;
for(int i=0;i<=m;i++)
for(int k=0;k<=m;k++)if(A.a[i][k])
for(int j=0;j<=m;j++)if(B.a[k][j])
res.a[i][j]=(res.a[i][j]+A.a[i][k]*B.a[k][j])%mod;
return res;
}
inline friend matrix Pow(matrix A,int b)
{
matrix res;
for(int i=0;i<=m;i++)res.a[i][i]=1;
for(;b;b>>=1,A=A*A)if(b&1)res=res*A;
return res;
}
}A,B;
inline int p(int i,int j){return j?i+a+b+1:i;}
int main()
{
//freopen("seq.in","r",stdin);
//freopen("seq.out","w",stdout);
n=getint(),a=getint(),b=getint(),mod=getint(),m=2*(a+b)+1;
for(int i=0;i<=a+b;i++)
{
c[i][0]=1;
for(int j=1;j<=i;j++)
c[i][j]=(c[i-1][j-1]+c[i-1][j])%mod;
}
A.a[0][p(0,0)]=1;
for(int i=0;i<=a+b;i++)
{
add(B.a[p(i,0)][p(i,0)],1),add(B.a[p(i,1)][p(i,0)],1);
for(int j=0;j<=i;j++)add(B.a[p(j,0)][p(i,1)],c[i][j]);
}
A=A*Pow(B,n);
for(int i=0,x=1;i<=a;i++,x=1ll*x*n%mod)
ans=(ans+c[a][i]*(((a-i)&1)?mod-1:1)%mod*x%mod*(A.a[0][p(a+b-i,0)]+A.a[0][p(a+b-i,1)])%mod)%mod;
cout<<ans<<'\n';
return 0;
}