图上定距离点对查找(邻接矩阵+矩阵快速幂+位运算优化)

yo 大家早上好、中午好、晚上好、凌晨好
欢迎来到本篇文章

简介

本文主要解决图上定距离点对查询的问题，此算法主要运用关系矩阵、矩阵快速幂、位运算，能以近 $O(n^2\log n)$ 的复杂度查找出所有存在距离为 $n$ 的路径的点对，并支持多次询问。

算法解释

关系矩阵

关系矩阵(matrix of a relation) 是对关系的一种刻画，即对于两个集合之间的某个关系，能清楚地表明此二集合的任意元素是否有此关系的数字矩阵。

关系矩阵是离散数学中的一个术语，其作用是为表示二元关系的一种方式，其外还有集合表达式和关系图可以用于表示二元关系。这里主要是利用了关系矩阵和关系图的关系，用邻接矩阵表示有向图（无向图可将边视为双向边）。在关系矩阵中用 $M_{i,j}=1$ 表示 $<i,j>$ 这一有序对存在，在关系图上表示为一条从 $i$ 连到 $j$ 的有向边，如二元关系：

$$\{<1,2>,<2,1>,<2,3>,<3,1>,<3,4>\}$$

可用关系图表示为：

而用关系矩阵则表示为：

$$\left\{\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 1& 0\\ 1& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right\}$$

没错，此时我们已经用只有0和1的邻接矩阵表示了一个有向图，而这是这个算法思想的基础。

矩阵快速幂

上面用了一个邻接矩阵表示了有向图上所有距离为1的点对，包括 $1\to 2,2\to 1,2\to 3,3\to 4$， 而下面我们便可以利用矩阵乘法将其转化成表示距离为 $n$ 的点对的矩阵

设 $R$ 与 $S$ 都为一个二元关系。$R$ 和 $S$ 的合成记作 $R\circ S$，其定义为 $R\circ S=\{<x,y>\mid \mathrm{\exists }t(<x,t>\in R\wedge <t,y>\in S)$，称为 $S$ 对 $R$ 的右复合。

离散数学中可用 $R^n$ 表示二元关系 $R$ 对自身的多次右复合，由右复合的定义可递推得出矩阵 $R$ 的 $n$ 次幂 $R^n$ 可表示图上距离为 $n$ 的点对集合。
则我们目标的矩阵 $M$ 可用矩阵快速幂由初始邻接矩阵，即关系矩阵 $R$ 算出。矩阵快速幂的代码这里就省略了，上面图中距离为6的点对最终算出矩阵 $M=R^6$ 为：

$$\left\{\begin{array}{cccc}2& 2& 1& 1\\ 3& 2& 2& 1\\ 2& 1& 1& 1\\ 0& 0& 0& 0\end{array}\right\}$$

位运算优化

一般的矩阵乘法大都要用到三层嵌套循环，其复杂度本身已经为 $O(n^3)$，相比之下还不如Floyd算法。
可观察到上面矩阵包含的元素其实不只0和1（此处每个元素其实表示点对间存在多少条距离为 $n$ 的路径），而我们需要确定的仅是确认是点对间是否存在定距离的路径，只需要0和1。矩阵乘法的定义有 $M_{i,j}=\sum _1^n{A}_{i,k}\times B_{k,j}$，而初始邻接矩阵中仅有0和1，则 $i$ 与 $j$ 间 $n$ 长度的路径不存在当且仅当 $R^{n-1}$ 中每对 $R_{i,k},R_{k,j}$ 均有0。
因此我们将原本用二维数组表示的邻接矩阵改为用一维数组表示。这里采用c++。

vector<unsigned LL> M(63, 0), MT(63, 0);
//MT用于存储矩阵M的转置，方便后续运算
//使用 unsigned long long 可使单个数表示的点达到64个。

这样便可以将矩阵乘法修改为矩阵位运算，整体代码如下：

/*
*日期：2023/3/30
*作者：AiHAn
*说明：此代码用于图上顶距离点对查找
*/

#include <bits/stdc++.h>
#define LL long long
#define pmm pair<vector<unsigned LL>, vector<unsigned LL> >

using namespace std;

struct edge
{
	int next, to;

}edges[10005];
int i = 1;
vector<int> head(55, 0); //链式前向星存图
vector<unsigned LL> MA(63, 0), MB(63, 0); //邻接矩阵
bitset<65> b;

void add(int u, int v)
{
	edges[i].to = v;
	edges[i].next = head[u];
	head[u] = i++;
}

//矩阵位运算
pmm matrixAnd(vector<unsigned LL> A, int size, vector<unsigned LL> AT)
{
	vector<unsigned LL> res1(55, 0), res2(55, 0);

	for (int j = 1; j <= size; j++)
		for (int k = 1; k <= size; k++)
			if (A[j] & AT[k]) //判断是否均不为0
			{
				res1[j] |= ((LL)1 << k); //1要用LL或unsigned LL表示，不然会变的不幸
				res2[k] |= ((LL)1 << j); //同时计算转置矩阵,方便下次运算
			}

	return make_pair(res1, res2);
}

//矩阵快速幂，和原版矩阵快速幂差不多，不加赘述
vector<unsigned LL> qpow(vector<unsigned LL> A, int size, vector<unsigned LL> AT, LL n)
{
	pmm res = make_pair(A, AT);
	vector<unsigned LL> res1(55, 0);
	for (int j = 1; j <= size; j++)
		res1[j] |= ((LL)1 << j);

	while (n)
	{
		if (n & 1)
			res1 = matrixAnd(res1, size, res.second).first;
		res = matrixAnd(res.first, size, res.second);
		n >>= 1;
	}

	return res1;
}

unsigned main()
{
	int n, m, u, v, l; //n为节点数，m为边数，l为目标距离
	cin >> n >> m >> l;
	for (int i = 0; i < m; i++)
	{
		cin >> u >> v;
		add(u, v);
		MA[u] |= ((LL)1 << v); //初始邻接矩阵
		MB[v] |= ((LL)1 << u); //初始邻接矩阵的转置
	}

	vector<unsigned LL> A = qpow(MA, n, MB, l);
	for (int j = 1; j <= n; j++) //输出
	{
		b = A[j];
		for (int k = 1; k <= n; k++)
			cout << b[k] << " ";
		cout << "\n";
	}
}

输入 #1：

4 5 3
1 2
2 1
2 3
3 1
3 4

输出 #1：

1 1 0 1 
1 1 1 0 
1 0 1 0 
0 0 0 0 

输入 #2：

4 5 6
1 2
2 1
2 3
3 1
3 4

输出：

1 1 1 1 
1 1 1 1 
1 1 1 1 
0 0 0 0 

高精度处理

上面仅使用单个 unsigned long long 数表示原邻接矩阵的一行，因此最多只能表示仅有64个节点的图，节点增多的话要写高精度。但由于涉及的运算仅有与运算和或运算，这里的高精度会比较容易写。代码如下：

vector<unsigned LL> operator & (vector<unsigned LL>& l, vector<unsigned LL>& r)
{
	vector<unsigned LL> res(10, 0);
	for (int i = 0; i < 10; i++)
		res[i] = l[i] & r[i];

	return res;
}

vector<unsigned LL> operator | (vector<unsigned LL>& l, vector<unsigned LL>& r)
{
	vector<unsigned LL> res(10, 0);
	for (int i = 0; i < 10; i++)
		res[i] = l[i] | r[i];

	return res;
}

而 num |= (LL)1 << k 等可改为 num[k % maxlength] = num[k % maxlength] | ((LL)1 << (k / maxlength)) 进行处理。

应用

luogu P1613

题意

给定一张边权均为 1 的有向图，图上距离为2的幂次的点对间可连一条长度也为 1 的边，求节点 1 到 $n$ 的最短路径长度。

解答

由于数据范围为 $n\le 50$，题解中大多都使用 Floyd + 倍增 这种 $O(n^3\log maxlongint)$ 的做法。这里我们用前文的算法来解决这道题，在进行快速幂时每次位运算后遍历一次矩阵，若 $M_{i,j}$ 为 1 就添加边，最后用 Dijstra 或 SPFA 等最短路算法求最短路即可 虽然这里作者用了最普通的BFS。最后时间复杂度仅为大概 $O(2n^2\log maxlongint)$。
下附AC代码。

代码

#include <bits/stdc++.h>
#define LL long long
#define pmm pair<vector<unsigned LL>, vector<unsigned LL> >

using namespace std;

struct edge
{
	int next, to;

}edges[10005]; //链式前向星存图
int i = 1, n;
vector<int> head(55, 0), vis(55, 1), cur(55, (int)1e9);
vector<vector<int> > evis(55, vector<int>(55, 1));
vector<unsigned LL> MA(55, 0), MAT(55, 0);
queue<int> q;
bitset<55> b;

void add(int u, int v)
{
	edges[i].to = v;
	edges[i].next = head[u];
	head[u] = i++;
}

//矩阵位运算
pmm matrixAnd(vector<unsigned LL> A, int size, vector<unsigned LL> BT)
{
	vector<unsigned LL> res1(55, 0), res2(55, 0);

	for (int j = 1; j <= size; j++)
		for (int k = 1; k <= size; k++)
			if (A[j] & BT[k])
			{
				res1[j] |= ((LL)1 << k);
				res2[k] |= ((LL)1 << j);
			}

	return make_pair(res1, res2);
}

//倍增
void qpow(vector<unsigned LL> A, int size, vector<unsigned LL> AT, LL N)
{
	pmm res = make_pair(A, AT);

	for(int j = 1; j <= N; j++)
	{
		res = matrixAnd(res.first, size, res.second);
		for (int k = 1; k <= size; k++) //每次乘方后遍历一次矩阵，添加新边
			for (int l = 1; l <= size; l++)
				if (evis[k][l] && (res.first[k] & ((LL)1 << l)))
				{
					add(k, l);
					evis[k][l] = 0;
				}
	}
}

void bfs(int p)
{
	if (p == n)
	{
		cout << cur[p];
		return;
	}

	for (int j = head[p]; j; j = edges[j].next)
		if (vis[edges[j].to])
		{
			cur[edges[j].to] = min(cur[edges[j].to], cur[p] + 1);
			q.push(edges[j].to);
			vis[edges[j].to] = 0;
		}

	q.pop();
	if (!q.empty())
		bfs(q.front());
}

unsigned main()
{
	int m, u, v;
	cin >> n >> m;
	for (int i = 0; i < m; i++)
	{
		cin >> u >> v;
		if (evis[u][v]) //去除重边降低BFS以时间复杂度
		{
			add(u, v);
			evis[u][v] = 0;		
			MA[u] |= ((LL)1 << v);
			MAT[v] |= ((LL)1 << u);
		}
	}

	qpow(MA, n, MAT, 31);

	q.push(1);
	vis[q.front()] = 0;
	cur[q.front()] = 0;
	bfs(q.front());
}

此算法纯属个人在离散课上忽然蹦出的想法
本文如有错误，还请大佬指正。