[NOI Online #3 提高组] 魔法值(矩阵快速幂+各种优化)

矩阵快速幂

相信其他巨佬讲矩阵快速幂的原理和过程都已经很详细了,这里我就着重讲一讲用裸的矩阵快速幂算法之后,如何优化。

优化1:优化矩阵乘法

这个好几位大佬讲过了,我也不讲了,主要就是把 O ( n 3 ) O(n^3) O(n3) 的矩阵乘法优化到 O ( n 2 ) O(n^2) O(n2)

优化2:倍增预处理矩阵

倍增预处理矩阵的核心思想就是:

2 k 2^k 2k 个转移矩阵相乘出来的矩阵预处理出来,并记录为 m k m_k mk,那么通过 m k = m k − 1 2 m_k=m_{k-1}^2 mk=mk12 就可以用 32 32 32 次矩阵乘法把每个 m i m_i mi 算出来。

那么对于一次询问为 x x x,也就是求 2 x 2^x 2x 个转移矩阵相乘出来的矩阵,我们就可以通过类似倍增一样的算法在 log ⁡ ( x ) log(x) log(x) 次计算后计算出答案。

优化3:bitset

我们可以发现,在计算矩阵乘法时,可以用 bitset优化。

比如,对于矩阵 A A A B B B,现在要计算 C = A × B C=A imes B C=A×B

那么根据定义, C i , j = xor ⁡ k = 1 n A i , k × B k , j C_{i,j}=operatorname{xor}_{k=1}^{n} A_{i,k} imes B_{k,j} Ci,j=xork=1nAi,k×Bk,j

形象地理解一下, C i , j C_{i,j} Ci,j 的值就是 A A A 矩阵的第 i i i 行和 B B B 矩阵的第 j j j 列一位一位地乘起来,再把这些乘出来的结果异或起来。

而且我们知道一个性质, A A A 矩阵和 B B B 矩阵中只可能出现 0 0 0 1 1 1

那么对于两个数 x , y ∈ { 0 , 1 } x,yin{0,1} x,y{0,1},必定有 x × y = x and ⁡ y x imes y=xoperatorname{and}y x×y=xandy

这时容易想到用 bitset维护矩阵乘法。就是把矩阵每一行所代表的 01 01 01 串存进 bitset<100>row[100]内,把矩阵的每一列所代表的 01 01 01 串存进 bitset<100>line[100]内。然后 A A A 矩阵的第 i i i 行和 B B B 矩阵的第 j j j 列一位一位地乘起来就是 A.row[i]&B.line[j],这样我们就能解决乘法的问题了。

但是有一个问题,我们得到了乘出来的序列后,怎么把它们一个一个异或起来呢?注意到这个序列是个 01 01 01 序列,那么如果这个序列中 1 1 1 的个数为奇数,它们异或起来就是 1 1 1,否则就是 0 0 0,所以异或的部分我们就可以维护了,即 C i , j = count ⁡ ( )   m o d   2 C_{i,j}=operatorname{count}()mod2 Ci,j=count()mod2

那么我们就顺利地把 O ( n 3 ) O(n^3) O(n3) 的矩阵乘法成功优化到 O ( n 3 64 ) O(frac{n^3}{64}) O(64n3)。(未加优化1)

优化4:询问离线

感觉最没用的优化。

显然,我们每次询问都算一次 pow() ⁡ operatorname{pow()} pow(),感觉会重复算了很多,所以不妨把询问给离线存下来,存进数组 q q q 里面,并对它进行排序。

这个排序可能会导致程序变慢而不是变快。

那么不妨假设我们现在已经知道了转移矩阵 t t t q i q_i qi 次方是多少,要求 q i + 1 q_{i+1} qi+1 的询问,也就是要求出 t q i + 1 t^{q_{i+1}} tqi+1。那么我们只用求出 t q i + 1 − q i t^{q_{i+1}-q_i} tqi+1qi,然后再乘上上次算出来的矩阵,就可以省下一些时间。

优化5:O2 优化和 IO 读写

最后贴一下代码:(没有加优化1卡过去的)

#include<bits/stdc++.h>

#define N 110
#define LV 32

using namespace std;

struct Matrix
{
	bitset<N>row[N];//bitset优化
	bitset<N>line[N];
	Matrix()
	{
		memset(row,0,sizeof(row));
		memset(line,0,sizeof(line));
	}
}turn[LV],st;

struct Query
{
	unsigned int x;
	int id;
}ask[N];

int n,m,q;
unsigned int f[N],ans[N];
bool tmp[N][N];

inline Matrix operator * (Matrix a,Matrix b)
{
	Matrix ans;
	bitset<N>now;
	for(register int i=1;i<=n;i++)
	{
		for(register int j=1;j<=n;j++)
		{
			now=a.row[i]&b.line[j];//上述的计算方法
			if(now.count()&1)
			{
				ans.row[i].set(j);
				ans.line[j].set(i);
			}
		}
	}
	return ans;
}

inline Matrix poww(unsigned int x)//倍增版的快速幂
{
	Matrix ans=st;
	for(int i=31;i>=0;i--)
	{
		if(x>=(1u<<i))
		{
			x-=(1u<<i);
			ans=turn[i]*ans;
		}
	}
	return ans;
}

inline void prework()//矩阵预处理
{
	for(register int i=1;i<LV;i++)
		turn[i]=turn[i-1]*turn[i-1];
}

inline bool cmp(const Query a,const Query b)
{
	return a.x<b.x;
}

inline unsigned int getans(const Matrix a)
{
	unsigned int ans=0;
	for(register int i=1;i<=n;i++)
		if(a.line[1][i])
			ans^=f[i];
	return ans;
}

inline int read()
{
	unsigned int x=0;
	char ch=getchar();
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1u)+(x<<3u)+(ch^'0');
		ch=getchar();
	}
	return x;
}

int main()
{
	n=read();m=read();q=read();
	for(register int i=1;i<=n;i++)
		st.row[i].set(i),st.line[i].set(i);
	for(register int i=1;i<=n;i++)
		f[i]=read();
	for(register int i=1;i<=m;i++)
	{
		const int u=read(),v=read();
		turn[0].row[u].set(v),turn[0].line[v].set(u);
		turn[0].row[v].set(u),turn[0].line[u].set(v);
	}
	prework();
   //下面是询问离线
	for(register int i=1;i<=q;i++)
	{
		ask[i].x=read();
		ask[i].id=i;
	}
	sort(ask+1,ask+q+1,cmp);
	Matrix last=poww(ask[1].x);
	ans[ask[1].id]=getans(last);
	for(register int i=2;i<=q;i++)
	{
		last=last*poww(ask[i].x-ask[i-1].x);
		ans[ask[i].id]=getans(last);
	}
	for(register int i=1;i<=q;i++)
		printf("%u
",ans[i]);
	return 0;
}
原文地址:https://www.cnblogs.com/ez-lcw/p/14448677.html