【矩阵快速幂倍增】P6569 [NOI Online #3 提高组] 魔法值|省选--优快云博客

本文链接：https://blog.youkuaiyun.com/he_zhidan/article/details/146116874

矩阵快速幂

P6569 [NOI Online #3 提高组] 魔法值

题目描述

H 国的交通由 $n$ 座城市与 $m$ 条道路构成，城市与道路都从 $1$ 开始编号，其中 $1$ 号城市是 H 国的首都。H 国中一条道路将把两个不同城市直接相连，且任意两个城市间至多有一条道路。

H 国是一个信奉魔法的国家，在第 $j$ 天， $i$ 号城市的魔法值为 $f_{i,j}$ 。H 国的魔法师已观测到第 0 天时所有城市的魔法值 $f_{i,0}$ ，且他们还发现，之后的每一天每个城市的魔法值，都将会变为所有与该城市直接相连的城市的前一天魔法值的异或值，即

$f_{x,j}=f_{v_1,j-1}\oplus f_{v_2,j-1}\oplus \cdots\oplus f_{v_k,j-1}$

其中 $j\ge 1$ ， $v_1,v_2,\cdots,v_k$ 是所有与 $x$ 号城市直接相连的城市， $\oplus$ 为异或运算。

现在 H 国的国王问了你 $q$ 个问题，对于第 $i$ （ $1\le i\le q$ ）个问题你需要回答：第 $a_i$ 天时首都的魔法值是多少。

输入格式

第一行三个用空格分隔的整数 $n, m, q$ ，表示城市数、道路数与问题数。

第二行 $n$ 个用空格分隔的整数，第 $i$ 个整数表示 $f_{i, 0}$ 。

接下来 $m$ 行，每行两个用空格分隔的正整数 $u, v$ ，表示一条连接 $u$ 号城市与 $v$ 号城市的道路。

接下来 $q$ 行每行一个整数，第 $i$ 行的整数表示 $a_i$ 。

输出格式

按顺序输出 $q$ 行每行一个整数，表示对应问题的答案。

输入输出样例 #1

输入 #1

输出 #1

说明/提示

数据规模与约定

对于 $20\%$ 的数据，满足 $a_i \leq 100$ 。
对于 $40\%$ 的数据，满足 $\leq 20$ 。
另有 $30\%$ 的数据，满足 $m=\frac{n(n-1)}{2}$ 。
对于 $100\%$ 的数据，满足 $\leq n,q \leq 100$ ， $\leq m \leq \frac{n(n-1)}{2}$ ， $1\leq a_i < 2^{32}$ ， $0\leq f_i < 2^{32}$ 。

说明

data provider ：@一扶苏一，数据有锅请联系她。如果被朴素的快速幂水过去了也请联系她。

异或和矩阵快速幂

普通矩阵乘法：各元素相乘后求和。本题各元素相乘后求异或和。本题也是遵守结合律的，故也可以用快速矩阵幂。
fi0就是pre,如果两个城市i,j向量则mat[i][j]=1,mat[j][i]=1 城市改成从0开始。
$pre \times mat^{ai}$ 便是答案。
时间复杂度：O(qnnnloga)
可以常数级优化：按a[i]从小到大处理（离线查询）,还是过不了。

倍增优化

预处理：
mats[i] = mat²ⁱ，预处理时间复杂度😮(nnnloga)
pre * mats[x] 注意：(1<<x) & a[i],任意查询，最多32个x。
pre一行N列，mats是N行N列，故一次矩阵运算 时间复杂度：O(nn)
一次查询时间复杂度：O(nnlogai)，总查询复杂度：O(qnnloga)

代码

核心代码

#include <iostream>
#include <sstream>
#include <vector>
#include<map>
#include<unordered_map>
#include<set>
#include<unordered_set>
#include<string>
#include<algorithm>
#include<functional>
#include<queue>
#include <stack>
#include<iomanip>
#include<numeric>
#include <math.h>
#include <climits>
#include<assert.h>
#include<cstring>
#include<list>
#include <array>

#include <bitset>
using namespace std;

template<class T1, class T2>
std::istream& operator >> (std::istream& in, pair<T1, T2>& pr) {
	in >> pr.first >> pr.second;
	return in;
}

template<class T1, class T2, class T3 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3>& t) {
	in >> get<0>(t) >> get<1>(t) >> get<2>(t);
	return in;
}

template<class T1, class T2, class T3, class T4 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3, T4>& t) {
	in >> get<0>(t) >> get<1>(t) >> get<2>(t) >> get<3>(t);
	return in;
}

template<class T = int>
vector<T> Read() {
	int n;
	scanf("%d", &n);
	vector<T> ret(n);
	for (int i = 0; i < n; i++) {
		cin >> ret[i];
	}
	return ret;
}

template<class T = int>
vector<T> Read(int n) {
	vector<T> ret(n);
	for (int i = 0; i < n; i++) {
		cin >> ret[i];
	}
	return ret;
}

template<int N = 1'000'000>
class COutBuff
{
public:
	COutBuff() {
		m_p = puffer;
	}
	template<class T>
	void write(T x) {
		int num[28], sp = 0;
		if (x < 0)
			*m_p++ = '-', x = -x;

		if (!x)
			*m_p++ = 48;

		while (x)
			num[++sp] = x % 10, x /= 10;

		while (sp)
			*m_p++ = num[sp--] + 48;
		AuotToFile();
	}
	void writestr(const char* sz) {
		strcpy(m_p, sz);
		m_p += strlen(sz);
		AuotToFile();
	}
	inline void write(char ch)
	{
		*m_p++ = ch;
		AuotToFile();
	}
	inline void ToFile() {
		fwrite(puffer, 1, m_p - puffer, stdout);
		m_p = puffer;
	}
	~COutBuff() {
		ToFile();
	}
private:
	inline void AuotToFile() {
		if (m_p - puffer > N - 100) {
			ToFile();
		}
	}
	char  puffer[N], * m_p;
};

template<int N = 1'000'000>
class CInBuff
{
public:
	inline CInBuff() {}
	inline CInBuff<N>& operator>>(char& ch) {
		FileToBuf();
		ch = *S++;
		return *this;
	}
	inline CInBuff<N>& operator>>(int& val) {
		FileToBuf();
		int x(0), f(0);
		while (!isdigit(*S))
			f |= (*S++ == '-');
		while (isdigit(*S))
			x = (x << 1) + (x << 3) + (*S++ ^ 48);
		val = f ? -x : x; S++;//忽略空格换行		
		return *this;
	}
	inline CInBuff& operator>>(long long& val) {
		FileToBuf();
		long long x(0); int f(0);
		while (!isdigit(*S))
			f |= (*S++ == '-');
		while (isdigit(*S))
			x = (x << 1) + (x << 3) + (*S++ ^ 48);
		val = f ? -x : x; S++;//忽略空格换行
		return *this;
	}
	template<class T1, class T2>
	inline CInBuff& operator>>(pair<T1, T2>& val) {
		*this >> val.first >> val.second;
		return *this;
	}
	template<class T1, class T2, class T3>
	inline CInBuff& operator>>(tuple<T1, T2, T3>& val) {
		*this >> get<0>(val) >> get<1>(val) >> get<2>(val);
		return *this;
	}
	template<class T1, class T2, class T3, class T4>
	inline CInBuff& operator>>(tuple<T1, T2, T3, T4>& val) {
		*this >> get<0>(val) >> get<1>(val) >> get<2>(val) >> get<3>(val);
		return *this;
	}
	template<class T = int>
	inline CInBuff& operator>>(vector<T>& val) {
		int n;
		*this >> n;
		val.resize(n);
		for (int i = 0; i < n; i++) {
			*this >> val[i];
		}
		return *this;
	}
	template<class T = int>
	vector<T> Read(int n) {
		vector<T> ret(n);
		for (int i = 0; i < n; i++) {
			*this >> ret[i];
		}
		return ret;
	}
	template<class T = int>
	vector<T> Read() {
		vector<T> ret;
		*this >> ret;
		return ret;
	}
private:
	inline void FileToBuf() {
		const int canRead = m_iWritePos - (S - buffer);
		if (canRead >= 100) { return; }
		if (m_bFinish) { return; }
		for (int i = 0; i < canRead; i++)
		{
			buffer[i] = S[i];//memcpy出错			
		}
		m_iWritePos = canRead;
		buffer[m_iWritePos] = 0;
		S = buffer;
		int readCnt = fread(buffer + m_iWritePos, 1, N - m_iWritePos, stdin);
		if (readCnt <= 0) { m_bFinish = true; return; }
		m_iWritePos += readCnt;
		buffer[m_iWritePos] = 0;
		S = buffer;
	}
	int m_iWritePos = 0; bool m_bFinish = false;
	char buffer[N + 10], * S = buffer;
};

class KMPEx
{
public:
	static vector<int> ZFunction(string s) {
		int n = (int)s.length();
		vector<int> z(n);
		z[0] = n;
		for (int i = 1, left = 0, r = 0; i < n; ++i) {
			if (i <= r) {//如果此if,r-i+1可能为负数
				z[i] = min(z[i - left], r - i + 1);
			}
			while ((i + z[i] < n) && (s[z[i]] == s[i + z[i]])) {
				z[i]++;
			}
			if (i + z[i] - 1 > r) left = i, r = i + z[i] - 1;
		}
		return z;//z[i] 表示S与其后缀S[i,n]的最长公共前缀(LCP)的长度
	}
	static int MinCyc(const string& str, int unit = 1) {
		const int N = str.length();
		auto z = ZFunction(str);
		auto Is = [&](int k) {
			for (int i = k; i < N; i <<= 1) {
				if (z[i] < min(N - i, i)) { return false; }
			}
			return true;
		};
		for (int k = unit; k < N; k += unit) { if (Is(k))return k; }
		return N;
	}

};

template<class T = long long, long long iDef = 0>
class CMatBase
{
public:
	CMatBase(const int N) :m_iN(N) {}
	// 矩阵乘法
	vector<vector<T>> multiply(const vector<vector<T>>& a, const vector<vector<T>>& b) {
		const int r = a.size(), c = b.front().size(), iK = a.front().size();
		assert(iK == b.size()); assert(iK == m_iN);
		vector<vector<T>> ret(r, vector<T>(c, iDef));
		for (int i = 0; i < r; i++)
		{
			for (int j = 0; j < c; j++)
			{
				for (int k = 0; k < iK; k++)
				{
					OnCal(ret[i][j], a[i][k], b[k][j]);
				}
			}
		}
		return ret;
	}
	// 矩阵快速幂
	vector<vector<T>> pow(const vector<vector<T>>& a, vector<vector<T>> b, T n) {
		vector<vector<T>> res = a;
		for (; n; n /= 2) {
			if (n % 2) {
				res = multiply(res, b);
			}
			b = multiply(b, b);
		}
		return res;
	}

	vector<vector<T>> CreateRow() {
		return vector<vector<T>>(1, vector<T>(m_iN));
	}
	vector<vector<T>> CreateUint() {
		vector<vector<T>> ret(m_iN, vector<T>(m_iN));
		for (int i = 0; i < m_iN; i++) { ret[i][i] = 1; }
		return ret;
	}
protected:
	virtual void OnCal(T& ans, const T& a1, const T& a2) = 0;
	const  int m_iN;
};
class CXorMat : public CMatBase<unsigned int, 0>
{
public:
	using CMatBase<unsigned int, 0>::CMatBase;
	virtual void OnCal(unsigned int& ans, const unsigned int& a1, const unsigned int& a2)
	{
		ans ^= a1 * a2;
	}
};
class Solution {
public:
	vector<unsigned int> Ans(vector<unsigned int>& f, vector<pair<int, int>> edge, vector< unsigned int >& que) {
		const int N = f.size();
		vector<vector<unsigned int>> pre(1, f);
		vector<vector<unsigned int>> mat(N, vector<unsigned int>(N));
		for (auto [i, j] : edge) {
			i--, j--;
			mat[i][j] = mat[j][i] = 1;
		}
		CXorMat matMul(N);
		vector < vector<vector<unsigned int>>> mats(32);
		mats[0] = mat;
		for (int i = 1; i < 32; i++) {
			mats[i] = matMul.multiply(mats[i - 1], mats[i - 1]);
		}
		vector<unsigned int> ans;
		for (const auto& q : que) {
			auto matAns = pre;
			for (int i = 0; i < 32; i++)
			{
				if ((1u << i) & q) {
					matAns = matMul.multiply(matAns, mats[i]);
				}
			}
			ans.emplace_back(matAns[0][0]);
		}
		return ans;
	}
};

int main() {
#ifdef _DEBUG
	freopen("a.in", "r", stdin);
#endif // DEBUG	
	ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
	int n, m, q;
	cin >> n >> m >> q;		
	auto f = Read<unsigned int>(n);
	auto edge = Read<pair<int, int>>(m);
	auto que = Read<unsigned int>(q);
#ifdef _DEBUG		
	/*printf("T=%d", T);
	Out(no1, ",no=");
	Out(strs, ",strs=");*/
	/*Out(edge, "edge=");
	Out(que, "que=");*/
#endif // DEBUG		
	auto res = Solution().Ans(f,edge,que);
	for (const auto& i : res) {
		cout << i << "\n";
	}
	return 0;
}

单元测试

vector<unsigned int> f, que;
		vector<pair<int, int>> edge;
		TEST_METHOD(TestMethod1)
		{
			f = { 0,0,1 }, edge = { {1,2},{1,3},{3,1} }, que = { 1 };
			auto res = Solution().Ans(f,edge,que);
			AssertEx({ 1 }, res);
		}

扩展阅读

我想对大家说的话
工作中遇到的问题，可以按类别查阅鄙人的算法文章，请点击《算法与数据汇总》。
学习算法：按章节学习《喜缺全书算法册》，大量的题目和测试用例，打包下载。重视操作
有效学习：明确的目标及时的反馈拉伸区（难度合适）专注
闻缺陷则喜(喜缺)是一个美好的愿望，早发现问题，早修改问题，给老板节约钱。
子墨子言之：事无终始，无务多业。也就是我们常说的专业的人做专业的事。
如果程序是一条龙，那算法就是他的是睛
失败+反思=成功成功+反思=成功

视频课程

先学简单的课程，请移步优快云学院，听白银讲师（也就是鄙人）的讲解。
https://edu.youkuaiyun.com/course/detail/38771
如何你想快速形成战斗了，为老板分忧，请学习C#入职培训、C++入职培训等课程
https://edu.youkuaiyun.com/lecturer/6176