淺談循環矩陣的乘法和快速冪


前言

前置技能:矩陣乘法,矩陣快速冪

當然你不會的話也不會點進來(滑稽)

今天上午的$HNOI$模擬賽中,$T1$是這么一道題目:

有一個長度為$n$的環,執行$s$次操作,在一次操作中,

對於每一個數,它變為它左邊的數乘上$l$以及它本身以及它右邊的數乘上$r$的和。

求最后每一個位置上的數是多少。(計算時左邊和右邊的數都是上一次的數)

最后結果模上$10^x$,$l,r,x$都為給定的常數

$n\leq1000,s\leq2^{30}$

很容易想到用矩陣快速冪來維護,假設我們現在有$4$個數字需要變換,設$f[i][j]$表示當前已經變換了$i$次,第$j$位上的數字是多少,有遞推式:
$$
\begin{aligned}
&\begin{bmatrix}
f[i-1][1]&f[i-1][2]&f[i-1][3]&f[i-1][4]
\end{bmatrix}
\times
\begin{bmatrix}
1&l&0&r\
r&1&l&0\
0&r&1&l\
l&0&r&1
\end{bmatrix}
\=
&\begin{bmatrix}
f[i][1]&f[i][2]&f[i][3]&f[i][4]
\end{bmatrix}
\end{aligned}
$$
但是這樣子做是$O(n^3log_2s)$,妥妥的超時(在沒有$O(wys)$優化的情況下)

這時我們發現,轉置矩陣是一個循環矩陣。

循環矩陣是什么?

下面內容引自百度百科,侵刪

​ 在線性代數中,循環矩陣是一種特殊形式的$Toeplitz$矩陣,它的行向量的每個元素都是前一個行向量各元素依次右移一個位置得到的結果。

簡單地說,就是以一行不斷向前/后變換的形式出現。循環矩陣有一些性質

假設矩陣$a,b$都是循環矩陣,那么:

  1. $a+b$是一個循環矩陣
  2. $a\times b$是一個循環矩陣

循環矩陣的矩陣快速冪

根據性質$2$,那么這個轉置矩陣無論自乘多少次,它還是一個循環矩陣,所以說,我們理論上只要知道第一行的最終形態,就可以知道整個矩陣了。那如何求出它的下一個形態呢?(轉置矩陣為$T$)

假設$g[i]$為當前矩陣第一行的第$i$個,$g'[i]$為下一個形態的矩陣的第$i$個。則有:
$$
\begin{aligned}
g'[1]=&T[1][1]\times T[1][1]+T[1][2]\times T[2][1]
+T[1][3]\times T[3][1]+T[1][4]\times T[4][1]
\
=&g[1]\times g[1]+g[2]\times g[4]+g[3]\times g[3]+g[4]\times g[2]
\end{aligned}
$$
這里可以得出這個結果是根據了循環矩陣的定義。

那么其他的求法也是雷同的,比如$g'[2]$,其他的交給讀者自己完成:
$$
\begin{aligned}
g'[2]=&T[1][1]\times T[1][2]+T[1][2]\times T[2][2]
+T[1][3]\times T[3][2]+T[1][4]\times T[4][2]
\
=&g[1]\times g[2]+g[2]\times g[1]+g[3]\times g[4]+g[4]\times g[3]
\end{aligned}
$$
於是我們得出了這樣的規律:
$$
g'[x]=\sum\limits_{(i+j-2)%n+1=x}g[i]g[j]
$$
很顯然,我們的初始矩陣$S$也是一個只有一行的循環矩陣,同樣適用於上面的規律,只不過是兩個不同的矩陣。
$$
tmp[x]=\sum\limits_{(i+j-2)%n+1=x}S[i]
T[j]
$$
其中$tmp$即本次矩陣乘法的結果矩陣

那么之前所講的例題再用矩陣快速冪就變成$O(n^2log_2s)$的了:

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using std::min; using std::max;
using std::swap; using std::sort;
typedef long long ll;

template<typename T>
void read(T &x) {
    int flag = 1; x = 0; char ch = getchar();
    while(ch < '0' || ch > '9') { if(ch == '-') flag = -flag; ch = getchar(); }
    while(ch >= '0' && ch <= '9') x = x * 10 + ch - '0', ch = getchar(); x *= flag;
}

const int N = 1e3 + 10;
int n, s, l, r, x, p = 1;
int S[N], T[N], tmp[N];

void mul(int S[], int T[]) {
	memset(tmp, 0, sizeof tmp);
	for(int i = 1; i <= n; ++i)
		for(int j = 1; j <= n; ++j)
			(tmp[(i + j - 2) % n + 1] += 1ll * S[i] * T[j] % p) %= p;
	memcpy(S, tmp, sizeof(tmp));
}//矩陣乘法

int main () {
	read(n), read(s), read(l), read(r), read(x);
	for(int i = 1; i <= x; ++i) p *= 10;
	for(int i = 1; i <= n; ++i) read(S[i]);
	T[1] = 1, T[2] = l, T[n] = r;
	for(; s; s >>= 1, mul(T, T)) if(s & 1) mul(S, T);//快速冪
	for(int i = 1; i <= n; ++i) printf("%d ", S[i]);
	return 0;
}


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM