xaybx^ay^bxayb十分不好维护,考虑怎么简化柿子
发现由于x+y=nx+y=nx+y=n,则为(n−y)ayb(n-y)^ay^b(n−y)ayb
用二项式定理暴力展开∑i=0aCia(−1)a−iniya+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=0aCia(−1)a−i\sum_{i=0}^aC_{i}^{a}(-1)^{a-i}∑i=0aCia(−1)a−i都是常量
那我们也就是要求出所有的ya+b−iy^{a+b-i}ya+b−i的和
考虑dpdpdp,令f[i][j][0/1]f[i][j][0/1]f[i][j][0/1]表示处理到第iii个数,结尾为0/1的111的数量的jjj次方和
则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=0jyk=∑k=0jCjkf[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=0jCkj10)\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)3logn)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;
}