【BZOJ5298】【CQOI2018】交错序列（矩阵快速幂优化dp）

传送门

$x^a{y}^b$十分不好维护，考虑怎么简化柿子
发现由于$x+y=n$，则为$(n-y{)}^a{y}^b$

用二项式定理暴力展开$$\sum _{i=0}^a{C}_i^a(-1{)}^{a-i}{n}^i{y}^{a+b-i}$$

发现对于所有柿子来说${\sum }_{i=0}^a{C}_i^a(-1{)}^{a-i}$都是常量
那我们也就是要求出所有的$y^{a+b-i}$的和

考虑$dp$，令$f[i][j][0/1]$表示处理到第$i$个数,结尾为0/1的$1$的数量的$j$次方和

则$f[i][j][0]=f[i-1][j][0]+f[i-1][j][1]$
$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][0],f[i][j][1])$$和转移矩阵

$$\left(\begin{array}{cc}1& {\sum }_{k=0}^j{C}_k^j\\ 1& 0\end{array}\right)$$

发现这样得出来就是$$(f[i+1][j][0],f[i+1][j][1])$$

矩阵快速幂优化即可
复杂度$O((a+b{)}^{3}logn)$

然后就可以拿到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;
}