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

传送门

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 (ny)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=0aCia(1)ainiya+bi

发现对于所有柿子来说 ∑ i = 0 a C i a ( − 1 ) a − i \sum_{i=0}^aC_{i}^{a}(-1)^{a-i} i=0aCia(1)ai都是常量
那我们也就是要求出所有的 y a + b − i y^{a+b-i} ya+bi的和

考虑 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[i1][j][0]+f[i1][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[i1][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) (11k=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;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值