矩阵乘法---快速幂

一、快速幂

给定三个正整数\(a,b,p\)，求出\(a^b%p\)的值。
代码如下：

int power(int a, int b, int p)
{
	int ans = 1 % p;
	while(b)
	{
		if(b & 1) ans = (long long) ans * a % p;
		b >>= 1;
		a = (long long) (a * a) % p;
	}
	return ans;
}

二、矩阵乘法+快速幂

给一个矩阵\(A\)，求出\(A^k\)的值（矩阵中所有元素都对\(mod\)取模）。
代码如下：

struct matrix
{
	int l, w;
	long long mat[111][111];
	matrix()
	{
		memset(mat, 0, sizeof mat);
	}//初始化
} A;
matrix create(int l, int w)
{
	matrix c;
	c.l = l, c.w = w;
	for(int i = 0; i < l; ++ i)
		for(int j = 0; j < w; ++ j)
			c.mat[i][j] = 1;
	return c;
}
matrix operator *(matrix a, matrix b)
{
	matrix c;
	c.l = a.l;
	c.w = b.w;
	for(int i = 0; i < a.l; ++ i)
		for(int j = 0; j < c.w; ++ j)
			for(int k = 0; k < a.w; ++ k)
				c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
	return c;
}
//calculate a ^ b % mod
matrix mat_power(matrix a, long long b) 
{
	matrix ans;
	ans.l = n, ans.w = n;
	if(b == 0) return create(n, n);
	for(int i = 0; i < n; ++ i) ans.mat[i][i] = 1; 
	while(b)
	{
		if(b & 1) ans = ans * a;//跟普通快速幂一致
		b >>= 1;
		a = a * a;
	}
	return ans;
}

例题：矩阵加速

网址：https://www.luogu.com.cn/problem/P1939

题目描述

已知一个数列\(a\)，它满足：
\(a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases}\)
求\(a\)数列的第\(n\)项对\(10^9+7\)取余的值。

输入格式

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

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

输出格式

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

输入输出样例

输入 #1

3
6
8
10

输出 #1

4
9
19

说明/提示

对于\(30%\)的数据\(n≤100\)；
对于\(60%\)的数据\(n≤2×10^7\)；
对于\(100%\)的数据\(1≤T≤100\)，\(1≤n≤2×10^9\)

---

构造矩阵Q即可。
代码如下：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const int mod = 1000000007;
int n;
struct matrix
{
	int l, w;
	long long mat[110][110];
	matrix()
	{
		memset(mat, 0, sizeof(mat));
	}
} f, Q;
matrix create(int l)
{
	matrix x;
	x.l = l; 
	x.w = l;
	for(int i = 0; i < l; ++ i) x.mat[i][i] = 1;
	return x;
}

matrix operator *(matrix a, matrix b)
{
	matrix c;
	c.l = a.l, c.w = b.w;
	for(int i = 0; i < c.l; ++ i)
		for(int j = 0; j < c.w; ++ j)
			for(int k = 0; k < a.w; ++ k)
				c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
	return c;
}
struct matrix mat_power(int b)
{
	matrix ans = f;
	while(b)
	{
		if(b & 1) ans = ans * Q;
		b >>= 1;
		Q = Q * Q;
	}
	return ans;
}
int main()
{
	f.l = 1, f.w = 3;
	f.mat[0][0] = 1, f.mat[0][1] = 1, f.mat[0][2] = 1;
	int T;
	scanf("%d", &T);
	while(T --)
	{
		Q.l = Q.w = 3;
		Q.mat[0][0] = Q.mat[0][1] = 0, Q.mat[0][2] = Q.mat[1][0] = 1, Q.mat[1][1] = Q.mat[1][2] = Q.mat[2][0] = 0, Q.mat[2][1] = Q.mat[2][2] = 1;
		scanf("%d", &n);
		if(n <= 3) puts("1");
		else printf("%d\n", mat_power(n - 3).mat[0][2]);
	}
	return 0;
}