【题解】矩阵加速

【题解】矩阵加速

小鱼干 167 2022-11-28

题目描述

已知一个数列 a，它满足：

ax={1x∈{1,2,3}ax−1+ax−3x≥4a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases}

求 a 数列的第 n 项对 109+710^9+7 取余的值。

输入格式

第一行一个整数 T，表示询问个数。

以下 T 行，每行一个正整数 n。

输出格式

每行输出一个非负整数表示答案。

输入输出样例

输入 #1

3
6
8
10

输出 #1

4
9
19

说明/提示

  • 对于 30% 的数据 n≤100n \leq 100；
  • 对于 60% 的数据 n≤2×107n \leq2 \times 10^7；
  • 对于 100% 的数据 1≤T≤100，1≤n≤2×1091 \leq T \leq 100，1 \leq n \leq 2 \times 10^9。

题目分析

题意很简单，求数列第n项的内容。但是注意n的范围，1≤n≤2×1091 \leq n \leq 2 \times 10^9 。即便是O(n)O(n)的递推，依然会超时。

从矩阵的角度观察递推式：f[n]=f[n-1]+f[n-3]

设矩阵A为：$ \begin{bmatrix}
a_{n-3} & a_{n-2} & a_{n-1}
\end{bmatrix} $。

[an−3an−2an−1]×[b1,1b1,2b1,3b2,1b2,2b2,3b3,1b3,2b3,3]=[an−2an−1an]\begin{bmatrix} a_{n-3} & a_{n-2} & a_{n-1} \end{bmatrix} \times \begin{bmatrix} b_{1,1} & b_{1,2} & b_{1,3} \\ b_{2,1} & b_{2,2} & b_{2,3} \\ b_{3,1} & b_{3,2} & b_{3,3} \end{bmatrix} = \begin{bmatrix} a_{n-2} & a_{n-1} & a_{n} \end{bmatrix}

那么意味着

an−3×b1,1+an−2×b2,1+an−1×b3,1=an−2an−3×b1,2+an−2×b2,2+an−1×b3,2=an−1an−3×b1,3+an−2×b2,3+an−1×b3,3=ana_{n-3}\times b_{1,1}+a_{n-2}\times b_{2,1}+ a_{n-1}\times b_{3,1}=a_{n-2}\\ a_{n-3}\times b_{1,2}+a_{n-2}\times b_{2,2}+ a_{n-1}\times b_{3,2}=a_{n-1}\\ a_{n-3}\times b_{1,3}+a_{n-2}\times b_{2,3}+ a_{n-1}\times b_{3,3}=a_{n}

那么加速矩阵B即可推出为

[001100011]\begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \end{bmatrix}

可发现，每乘一次矩阵B，即可向后推一项。那么求第n项就变成了求 A×BnA\times B^n 。我们可以采用矩阵快速幂的方式在O(logn)O(logn)的时间复杂度内求解。

代码实现

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const int N=5;
const int M=1e9+7;
struct node{
	ll a[N][N]={0};
	int row,col;
};
node I;//单位矩阵

node matrixMins(node a,node b){//矩阵乘法
	node c;//答案矩阵
	c.row=a.row;
	c.col=b.col;
	int n=c.row,p=c.col,m=a.col;
	//计算矩阵乘法
	for(int i=1;i<=n;i++){
		for(int j=1;j<=p;j++){
			for(int k=1;k<=m;k++){
				c.a[i][j]+=a.a[i][k]*b.a[k][j]%M;
				c.a[i][j]%=M;
			}
		}
	}
	return c;
}
node matrixPow(node a,ll k){//矩阵的幂次方
	if(k==0){// 0次方
		return I;//矩阵的0次方是单位矩阵
	}
	node t=matrixPow(a,k/2);//求 a^{n/2} 次方
	if(k&1){//判断k是否是奇数
		return matrixMins(matrixMins(t,t),a);
	}else{//k是偶数
		return matrixMins(t,t);
	}
}
int main(){
	int t;
	node a;
	ll n;
	cin>>t;

	//处理 加速递推矩阵
	a=node{
		{
			{0},
			{0,0,0,1},
			{0,1,0,0},
			{0,0,1,1}
		},3,3
	};
	//处理 单位矩阵
	I=node{
		{
			{0},
			{0,1,0,0},
			{0,0,1,0},
			{0,0,0,1}
		},3,3
	};
	//处理数列初始值 [1 1 2 3]
	node tt={
		{
			{0},
			{0,0,1,1}
		},1,3
	};

	while(t--){
		cin>>n;
		node tmp=matrixPow(a,n);//计算A^n
		node ans=matrixMins(tt,tmp);
		cout<<ans.a[1][1]<<endl;		
	}
	return 0;
}
算法5 矩阵快速幂 矩阵