[SDOI2008]递归数列 题解

矩乘板子题

(k le 15)(a_{i}= { extstyle sum_{j=i}^{1}}c_{j}*a_{i-j} (i>k)) (这一看就很矩乘) 考虑矩阵加速。

题目让求一段区间的和,可以转化为两前缀和相减的形式,同时,让矩阵边长加上 (1) , 多加的 (1) 用来储存前缀和

状态矩阵 (f_{i}) 存储:

[egin{bmatrix} a_{i}&a_{i+1} &a_{i+2} &··· &a_{i+k} &S_{i-1} end{bmatrix}]

我们想让他转移到 (f_{i+1}) 即为:

[egin{bmatrix} a_{i+1}&a_{i+2} &a_{i+3} &··· &a_{i+k+1} &S_{i} end{bmatrix}]

很容易得到,转移矩阵 (A) 为:

[egin{bmatrix} 0 &0 &0 &··· &c_{k} &1 \ 1 &0 &0 &··· &c_{k-1} &0 \ 0 &1 &0 &··· &c_{k-2} &0 \ ··· &··· &··· &··· &··· &··· \ 0 &0 &0 &··· &c_{1} &0 \ 0 &0 &0 &··· &0 &1 end{bmatrix}]

那矩乘的柿子就有了:

[egin{bmatrix} a_{i}&a_{i+1} &a_{i+2} &··· &a_{i+k} &S_{i-1} end{bmatrix} egin{bmatrix} 0 &0 &0 &··· &c_{k} &1 \ 1 &0 &0 &··· &c_{k-1} &0 \ 0 &1 &0 &··· &c_{k-2} &0 \ ··· &··· &··· &··· &··· &··· \ 0 &0 &0 &··· &c_{1} &0 \ 0 &0 &0 &··· &0 &1 end{bmatrix} = egin{bmatrix} a_{i+1}&a_{i+2} &a_{i+3} &··· &a_{i+k+1} &S_{i} end{bmatrix} ]

举个栗子,对于 (k=4) 时,有:

[egin{bmatrix} a_{i}& a_{i+1}& a_{i+2}& a_{i+3}& S_{i-1} end{bmatrix} egin{bmatrix} 0& 0& 0& c_{4}& 1\ 1& 0& 0& c_{3}& 0\ 0& 1& 0& c_{2}& 0\ 0& 0& 1& c_{1}& 0\ 0& 0& 0& 0&1 end{bmatrix} = egin{bmatrix} a_{i+1}& a_{i+2}& a_{i+3}& a_{i+4}& S_{i} end{bmatrix}]

(S_{0}) 计算到 (S_{n}) 需要进行矩阵乘法需要使 (f_{0})(n)(A),

及要计算 (f_{0}*A^{n}) 计算 (A^{n}) 用矩阵快速幂加速运算即可。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 17
#define LL long long 
using namespace std;

LL k,b[N],c[N],n,m,mod;

struct matrix
{
	LL a[N][N];
	void clear()//矩阵清空
	{
		memset(a,0,sizeof(a));
	}
	void build()//建立单位矩阵
	{
		memset(a,0,sizeof(a));
		for(register int i=1;i<=k;i++)
			a[i][i]=1;
	}
	void mat_mod()//矩阵取模
	{
		for(register int i=1;i<=k;i++)
			for(register int j=1;j<=k;j++)
					a[i][j]=a[i][j]%mod;
	}
	void print()//矩阵输出
	{
		for(register int i=1;i<=k;i++)
		{
			for(register int j=1;j<=k;j++)
				printf("%lld ",a[i][j]);
			printf("
");
		}
	}
};

matrix operator *(const matrix x,const matrix y)//矩阵乘法
{
	matrix ans;
	ans.clear();
	for(register int i=1;i<=k;i++)
		for(register int j=1;j<=k;j++)
			for(register int p=1;p<=k;p++)
				ans.a[i][j]=(ans.a[i][j]+x.a[i][p]*y.a[p][j]%mod)%mod;
	return ans;
}

inline LL qr()
{
	char ch=0;LL w=1,x=0;
	while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
	return x*w;
}

matrix poww(matrix x,LL opl)//矩阵快速幂
{
	matrix ans;
	ans.build();//创建单位矩阵
	while(opl)
	{
		if(opl&1)
			ans=ans*x;
		x=x*x;
		opl>>=1;
	}
	return ans;
}

int main()
{
	k=qr();
	for(register int i=1;i<=k;i++)
		b[i]=qr();
	for(register int i=1;i<=k;i++)
		c[i]=qr();
	n=qr();m=qr();mod=qr();
	k++;

	matrix A;//A为转化矩阵
	A.clear();
	for(register int i=2;i<k;i++)//给A填数
		A.a[i][i-1]=1;
	for(register int i=1;i<k;i++)
		A.a[i][k-1]=c[k-i];
	A.a[1][k]=1;
	A.a[k][k]=1;
	A.mat_mod();

	matrix f1,f2;
	f1.clear();//f1,f2为状态矩阵
	f2.clear();
	for(register int i=1;i<k;i++)//给f1,f2填数
		f1.a[1][i]=f2.a[1][i]=b[i];
	f1.mat_mod();
	f2.mat_mod();

	matrix jc1=poww(A,n-1ll);
	matrix jc2=poww(A,m);//A^m
	f1=f1*jc1;//计算前缀和
	f2=f2*jc2;

	printf("%lld
",((f2.a[1][k]-f1.a[1][k])%mod+mod)%mod);//前缀和减一下即为最终答案
	return 0;
}
原文地址:https://www.cnblogs.com/isonder/p/14375830.html