矩阵快速幂加速递推

P3390 【模板】矩阵快速幂#

矩阵乘法：

$c_{i, j} = \sum a_{i, k} * b_{k, j}$

矩阵快速幂：与普通快速幂类似。

特殊地，定义 $A^{0}$ 为单位矩阵 $I = [\begin{matrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 1 \end{matrix}]$ 。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;
using i64 = long long;

const int N = 110, mod = 1e9 + 7;

struct matrix {
	matrix() { memset(a, 0, sizeof(a)); }
	int a[N][N];
};

int n;
i64 k;

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	for (int k = 1; k <= n; k++)
		for (int i = 1; i <= n; i++)
			for (int j = 1; j <= n; j++)
				c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
	return c;
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	cin >> n >> k;
	
	matrix ans, b;
	
	for (int i = 1; i <= n; i++) ans.a[i][i] = 1;
	
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++)
			cin >> b.a[i][j];
	
	while (k) {
		if (k & 1) ans = ans * b;
		b = b * b;
		k >>= 1;
	}
	
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) cout << ans.a[i][j] << ' ' ;
		cout << '\n';
	}
	return 0;
}

POJ-3070/P1962 斐波那契数列#

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;
using i64 = long long;

const int N = 3, mod = 1e9 + 7;

struct matrix {
	matrix() { memset(a, 0, sizeof(a)); }
	int a[N][N];
};

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	
	for (int i = 1; i < N; i++)
		for (int j = 1; j < N; j++)
			for (int k = 1; k < N; k++)
				c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
	
	return c;
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	matrix ans, b;
	i64 n;
	
	cin >> n;
	
	for (int i = 0; i < N; i++) ans.a[i][i] = 1;
	
	b.a[1][1] = 1;
	b.a[1][2] = 1;
	b.a[2][1] = 1;
	b.a[2][2] = 0;
	
	while (n) {
		if (n & 1) ans = ans * b;
		b = b * b;
		n >>= 1;
	}
	
	cout << ans.a[1][2] << '\n';
	return 0;
}

P1939 【模板】矩阵加速（数列）#

本题是斐波那契数列的进阶版。

因为 $a_{x}$ 的计算涉及三项： $a_{x - 1}$ ， $a_{x - 2}$ ， $a_{x - 3}$ 。

所以我们有

$[\begin{matrix} a_{x} & a_{x - 1} & a_{x - 2} \end{matrix}] = [\begin{matrix} a_{x - 1} & a_{x - 2} & a_{x - 3} \end{matrix}] \times A$

$A$ 应为：

$A = [\begin{matrix} a & b & c \\ d & e & f \\ g & h & i \end{matrix}]$

按照套路，可以解得：

$A = [\begin{matrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix}]$

所以原问题转化为 $[\begin{matrix} a_{x + 2} & a_{x + 1} & a_{x} \end{matrix}] = [\begin{matrix} a_{3} & a_{2} & a_{1} \end{matrix}] \times A^{x}$

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;

const int N = 4, mod = 1e9 + 7; 

struct matrix {
	matrix() { memset(a, 0, sizeof(a));	}
	int a[N][N];
};

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	
	for (int i = 1; i < N; i++)
		for (int j = 1; j < N; j++)
			for (int k = 1; k < N; k++)
				c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
				
	return c;
}

void solve() {
	int n;
	cin >> n;
	
	matrix ans, b;
	
	for (int i = 0; i < N; i++) ans.a[i][i] = 1;
	
	b.a[1][1] = 1;
	b.a[1][2] = 1;
	b.a[2][3] = 1;
	b.a[3][1] = 1;
	
	while (n) {
		if (n & 1) ans = ans * b;
		b = b * b;
		n >>= 1;
	}
	
	memset(b.a, 0, sizeof(b.a));
	b.a[1][1] = 1;
	b.a[1][2] = 1;
	b.a[1][3] = 0;
	
	ans = b * ans;
	
	cout << ans.a[1][3] << '\n';
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	int T;
	
	cin >> T;
	while (T--) solve();
	
	return 0;
}

POJ 3233/UVA 11149 Power of Matrix#

以UVA的题目为例~~，POJ上的题目就是把solve()函数执行一遍即可。~~：

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;

const int N = 90, mod = 10;	

struct matrix {
	matrix() { memset(a, 0, sizeof(a));	}
	int a[N][N];
};

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	
	for (int i = 1; i < N; i++)
		for (int j = 1; j < N; j++)
			for (int k = 1; k < N; k++)
				c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j]) % mod;
	
	return c;
}

int n, k;

void solve() {
	matrix ans, b;
	
	for (int i = 1; i <= (n << 1); i++) ans.a[i][i] = 1;
	
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++) {
			int x;
			cin >> x;
			x %= 10;
			b.a[i][j] = x;
			b.a[i][j + n] = 0;
			b.a[i + n][j] = x;
			if (i == j) b.a[i + n][j + n] = 1;
			else b.a[i + n][j + n] = 0;
		}
	}
	
	while (k) {
		if (k & 1) ans = ans * b;
		b = b * b;
		k >>= 1;
	}
	
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j < n; j++) {
			cout << ans.a[i + n][j] << ' ';
		}
		cout << ans.a[i + n][n] << '\n';
	}
	
	cout << '\n';
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	while (cin >> n >> k, (n && k)) solve();
	return 0;
}

P4910 帕秋莉的手环#

~~吐槽一句：题面写得太好了！！！~~

建议看翻译。

本题先使用普通DP进行求解：

设 $f [i] [0 / 1]$ 为前 $i - 1$ 位都已经选择好，且第 $i$ 位染上了颜色 $0 / 1$ 的方案数。（设 $1$ 为黑色， $0$ 为白色。

于是有了：

$f [i] [1] = f [i - 1] [0]$ ,

$f [i] [0] = f [i - 1] [0] + f [i - 1] [1]$

而第一个珠子分为两种情况，黑的或白的，分别dp就可以了。

这种暴力dp可以拿 48 分~~，dalao们有的可以 AC~~。

我们通过运算，不难发现，答案为：

设 $f i b [i]$ 为斐波那契数列的第 $i$ 项。

$a n s = f i b [n - 1] \times 2 + f i b [n]$ 。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;
using i64 = long long;

const int N = 3, mod = 1e9 + 7;

struct matrix {
	matrix() { memset(a, 0, sizeof(a));	}
	int a[N][N];
};

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	for (int i = 1; i < N; i++)
		for (int j = 1; j < N; j++)
			for (int k = 1; k < N; k++)
				c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
	return c;
}

void solve() {
	i64 n;
	cin >> n;
	n--;
	
	matrix a, res;
	
	for (int i = 1; i < N; i++) res.a[i][i] = 1;
	
	a.a[1][1] = 1;
	a.a[1][2] = 1;
	a.a[2][1] = 1;
	a.a[2][2] = 0;
	
	while (n) {
		if (n & 1) res = res * a;
		a = a * a;
		n >>= 1;
	}
	cout << (res.a[1][2] * 2 % mod + res.a[1][1]) % mod << '\n';
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	int T;
	cin >> T;
	while (T--) solve(); 
	return 0;
}

P2886 [USACO07NOV]Cow Relays G#

题意：
给定一张 $T$ 条边的无向连通图，求从 $S$ 到 $E$ 经过 $N$ 条边的最短路长度。

用邻接矩阵 $M$ 表示一个图， $M (i, j)$ 是其第 $i$ 行第 $j$ 列元素，表示点 $i$ 和点 $j$ 的直接关系，若点 $i$ 和点 $j$ 直连，令 $M [i, j]$ 等于 $i j$ 的权值，否则等于无穷大（ $\infty$ ）；当 $i$ 等于 $j$ 时令 $M [i, i] = \infty$ 。

定义矩阵的广义乘法：

$M \times M = {min}_{a = 1}^{N} [M (i, a) + M (a, j)]$ ，

也就是把普通的矩阵乘法从求和改成了取最小值，把内部项相乘改为了相加。

定理：计算邻接矩阵的广义幂 $G = M^{n}$ ， $G (i, j)$ 的值是从 $i$ 到 $j$ 经过 $n$ 条边（或 $n - 1$ 个点）的最短路径长度。

以上摘自《算法竞赛》。

注意：在做矩阵乘法时，一定要处理好上界(代码中为 $i d x$ )。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

using namespace std;

const int N = 1010, M = 110;

int n, m;
int S, T;

int ver[N], idx;

int add(int a) {
	if (ver[a]) return ver[a];
	ver[a] = ++idx;
	return ver[a]; 
}

struct matrix {
	matrix() { memset(a, 0x3f, sizeof(a));	}
	int a[M][M];
};

matrix operator*(matrix& a, matrix& b) {
	matrix c;
	for (int k = 1; k <= idx; k++)
		for (int i = 1; i <= idx; i++)
			for (int j = 1; j <= idx; j++)
				c.a[i][j] = min(c.a[i][j], a.a[i][k] + b.a[k][j]);
	return c;
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	
	matrix ans, g;
	cin >> n >> m >> S >> T;
	
	for (int i = 1; i <= m; i++)
	{
		int w, a, b, t1, t2;
		cin >> w >> a >> b;
		t1 = add(a);
		t2 = add(b);
		g.a[t1][t2] = g.a[t2][t1] = w;
	}
	
	n--;
	ans = g;
	
	while (n) {	
		if (n & 1) ans = ans * g;
		g = g * g;
		n >>= 1;
	}
	
	cout << ans.a[add(S)][add(T)] << '\n';
	return 0;
}

[参考文献]

《算法竞赛》

posted @ 2022-12-03 00:17 SunnyYuan 阅读(10) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 线段树的基本操作与线段树合并

· Aho-Corasick自动机（AC自动机）

· 【模板】矩阵快速幂

· 矩阵快速幂笔记

· 数学基础-快速幂、快速乘、矩阵快速幂

公告

昵称： SunnyYuan
园龄： 2年7个月
粉丝： 15
关注： 5

+加关注

2025年3月

日

一

二

三

四

五

六

Loading

SunnyYuan的博客

矩阵快速幂加速递推

P3390 【模板】矩阵快速幂#

POJ-3070/P1962 斐波那契数列#

P1939 【模板】矩阵加速（数列）#

POJ 3233/UVA 11149 Power of Matrix#

P4910 帕秋莉的手环#

P2886 [USACO07NOV]Cow Relays G#

公告

搜索

随笔档案

阅读排行榜

推荐排行榜

最新评论