x
a
y
b
x^ay^b
xayb十分不好维护,考虑怎么简化柿子
发现由于
x
+
y
=
n
x+y=n
x+y=n,则为
(
n
−
y
)
a
y
b
(n-y)^ay^b
(n−y)ayb
用二项式定理暴力展开 ∑ i = 0 a C i a ( − 1 ) a − i n i y a + b − i \sum_{i=0}^aC_{i}^{a}(-1)^{a-i}n^iy^{a+b-i} i=0∑aCia(−1)a−iniya+b−i
发现对于所有柿子来说
∑
i
=
0
a
C
i
a
(
−
1
)
a
−
i
\sum_{i=0}^aC_{i}^{a}(-1)^{a-i}
∑i=0aCia(−1)a−i都是常量
那我们也就是要求出所有的
y
a
+
b
−
i
y^{a+b-i}
ya+b−i的和
考虑 d p dp dp,令 f [ i ] [ j ] [ 0 / 1 ] f[i][j][0/1] f[i][j][0/1]表示处理到第 i i i个数,结尾为0/1的 1 1 1的数量的 j j j次方和
则
f
[
i
]
[
j
]
[
0
]
=
f
[
i
−
1
]
[
j
]
[
0
]
+
f
[
i
−
1
]
[
j
]
[
1
]
f[i][j][0]=f[i-1][j][0]+f[i-1][j][1]
f[i][j][0]=f[i−1][j][0]+f[i−1][j][1]
f
[
i
]
[
j
]
[
1
]
=
(
y
+
1
)
j
=
∑
k
=
0
j
y
k
=
∑
k
=
0
j
C
j
k
f
[
i
−
1
]
[
k
]
[
0
]
f[i][j][1]=(y+1)^j=\sum_{k=0}^{j}y^k=\sum_{k=0}^{j}C_{j}^{k}f[i-1][k][0]
f[i][j][1]=(y+1)j=∑k=0jyk=∑k=0jCjkf[i−1][k][0]
发现很像矩阵乘法的形式
考虑构造答案矩阵 ( f [ i ] [ j ] [ 0 ] , f [ i ] [ j ] [ 1 ] ) (f[i][j][0],f[i][j][1]) (f[i][j][0],f[i][j][1])和转移矩阵
( 1 ∑ k = 0 j C k j 1 0 ) \left(\begin{array}{cc} 1 & \sum_{k=0}^{j}C_{k}^{j} \\ 1 & 0 \end{array}\right) (11∑k=0jCkj0)
发现这样得出来就是 ( f [ i + 1 ] [ j ] [ 0 ] , f [ i + 1 ] [ j ] [ 1 ] ) (f[i+1][j][0],f[i+1][j][1]) (f[i+1][j][0],f[i+1][j][1])
矩阵快速幂优化即可
复杂度
O
(
(
a
+
b
)
3
l
o
g
n
)
O((a+b)^3logn)
O((a+b)3logn)
然后就可以拿到40分的好成绩了
但发现会妥妥T掉
考虑怎么优化
1:发现转移矩阵有用的只有四个三角部分,优化一半常数(没写)
2:取模优化
3:发现m只有1e8,可以先爆乘,等过大时再取模(跑的贼快)
#include<bits/stdc++.h>
using namespace std;
#define ll long long
inline int read(){
char ch=getchar();
int res=0,f=1;
while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
while(isdigit(ch))res=res*10+(ch^48),ch=getchar();
return res*f;
}
const ll maxv=1LL<<62;
#define I64(x) static_cast<long long>(x)
int n,a,b,mod,hs,siz;
const int N=185;
inline int add(int a,int b){
return (a+b>=mod)?(a+b-mod):(a+b);
}
inline int mul(int a,int b){
return 1ll*a*b%mod;
}
struct matrix{
int n,m,a[N][N];
matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
matrix operator*(const matrix& rhs) const{
matrix res(n,rhs.m);
for(int i=0;i<n;++i)
for(int j=0;j<rhs.m;++j){
ll sum=0;
for(int k=0;k<m;++k){
sum+=I64(a[i][k])*rhs.a[k][j];
if(sum>=maxv)sum%=mod;
}
res.a[i][j]=sum%mod;
}
return res;
}
};
int c[N][N];
ll pw[N],anc;
int main(){
n=read(),a=read(),b=read(),mod=read();
hs=a+b+1,siz=hs*2;
c[0][0]=1;
for(int i=1;i<=hs;i++){
c[i][0]=1;
for(int j=1;j<=i;j++)
c[i][j]=add(c[i-1][j-1],c[i-1][j]);
}
pw[0]=1;
for(int i=1;i<=a;i++)pw[i]=pw[i-1]*n%mod;
matrix res(siz,siz);
for(int i=0;i<hs;i++){
res.a[i][i]=1,res.a[i+hs][i]=1;
for(int j=i;j<hs;j++)
res.a[i][j+hs]=c[j][i];
}
matrix ans(1,siz);
ans.a[0][0]=1;
for(;n;n>>=1,res=res*res){
if(n&1)ans=ans*res;
}
for(int i=0;i<=a;i++){
int p=a+b-i;
ll fc=c[a][i]*(((a-i)&1)?-1:1)*pw[i]%mod;
ll t=(ans.a[0][p]+ans.a[0][p+hs])%mod;
anc+=t*fc,anc%=mod;
}
cout<<anc;
}