矩阵快速幂
相信其他巨佬讲矩阵快速幂的原理和过程都已经很详细了,这里我就着重讲一讲用裸的矩阵快速幂算法之后,如何优化。
优化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=mk−12 就可以用 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\times 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}\times 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,y\in\{0,1\} x,y∈{0,1},必定有 x × y = x and y x\times y=x\operatorname{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}()\bmod2 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+1−qi,然后再乘上上次算出来的矩阵,就可以省下一些时间。
优化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\n",ans[i]);
return 0;
}