關於循環矩陣的一些理解和應用

問題引入

給定一個長度為\(n\)的環，每次對於一個位置上的數，加上它左邊的數乘上\(L\)，再加上右邊的數乘上\(R\)。\(L,R\)是給定的常數，問執行\(s\)次之后，這個環上的每個位置上的數是多少

（計算時左邊和右邊的數都是上一次的數）

如果\(n \leq 100, s \leq 2^{30}\)

可以想到矩陣快速冪

構造矩陣

\[\begin{aligned} &\begin{bmatrix}a_1 & a_2 & \cdots & a_n\end{bmatrix} \times \begin{bmatrix} 1 & l & 0 & 0 & \cdots & 0 & r \\ r & 1 & l & 0 & \cdots & 0 & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ l & 0 & \cdots & \cdots & \cdots & r & 1 \end{bmatrix} \\ =\;& \begin{bmatrix}a_1' & a_2' & \cdots & a_n'\end{bmatrix} \end{aligned} \]

於是開心地快速冪了

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define RG register
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define clear(x, y) memset(x, y, sizeof(x))

inline int read()
{
	int data = 0, w = 1; char ch = getchar();
	while(ch != '-' && (!isdigit(ch))) ch = getchar();
	if(ch == '-') w = -1, ch = getchar();
	while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
	return data * w;
}

const int maxn(1010);
int n, s, L, R, _x, Pow[10], a[maxn], Mod;

namespace $
{
	const int N(110);
	struct Matrix
	{
		int n, m, a[N][N];
		Matrix(int x, int y) : n(x), m(y) { memset(a, 0, sizeof a); }
		int *operator [] (const int &id) { return a[id]; }
		const int *operator [] (const int &id) const { return a[id]; }
	};

	Matrix operator * (const Matrix &a, const Matrix &b)
	{
		Matrix c(a.n, b.m);
		for(RG int i = 0; i < a.n; i++)
			for(RG int j = 0; j < a.m; j++)
				for(RG int k = 0; k < b.m; k++)
					c[i][k] = (1ll * c[i][k] + 1ll * a[i][j] * b[j][k]) % Mod;
		return c;
	}

	int main()
	{
		Matrix S(1, n), T(n, n);
		for(RG int i = 0; i < n; i++) S[0][i] = a[i];
		for(RG int i = 0; i < n; i++)
		{
			int pre = (i - 1 + n) % n, suc = (i + 1) % n;
			T[i][i] = 1, T[pre][i] = L, T[suc][i] = R;
		}
		while(s)
		{
			if(s & 1) S = S * T;
			T = T * T; s >>= 1;
		}
		for(RG int i = 0; i < n; i++) printf("%d ", S[0][i]);
		return 0;
	}
}

int main()
{
	Pow[0] = 1;
	for(RG int i = 1; i < 10; i++) Pow[i] = Pow[i - 1] * 10;
	n = read(), s = read(), L = read(), R = read(), _x = read();
	for(RG int i = 0; i < n; i++) a[i] = read();
	Mod = Pow[_x];
	return $::main();
}

解決

可是：

\(n \leq 1000, s \leq 2^{30}\)

\(\text{TLE}\)

但仔細看這個矩陣可以發現一個特點：它的每一行都是上一行往右移動一位得到的

不但如此，這個矩陣無論自乘多少次，都滿足這個性質，所以理論上只需要維護第一行就好了

令\(k[i]\)表示矩陣第一行的第\(i\)個，那么有（假設是一個\(4 * 4\)的矩陣

\[\begin{aligned} k[1]'&=a[1][1]\times a[1][1] + a[1][2]\times a[2][1] + a[1][3]\times a[3][1] + a[1][4]\times a[4][1] \\ &=k[1]\times k[1] + k[2]\times k[4] + k[3]\times k[3] + k[4]\times k[2] \end{aligned} \\ \therefore k[2] = k[1]\times k[2] + k[2]\times k[1] + k[3]\times k[4] + k[4]\times k[3] \]

可以看出一些規律：

\[k[x]=\sum_{(i+j-2) \% n + 1 = x} k[i] \times k[j] \]

這樣復雜度就是\(O(n^2log\;s)\)了

解決啦！！！

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define RG register
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define clear(x, y) memset(x, y, sizeof(x))

inline int read()
{
	int data = 0, w = 1; char ch = getchar();
	while(ch != '-' && (!isdigit(ch))) ch = getchar();
	if(ch == '-') w = -1, ch = getchar();
	while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
	return data * w;
}

const int maxn(1010);
int n, s, L, R, _x, Mod = 1;
int S[maxn], T[maxn];

void Mul(int S[], int T[])
{
	static int c[maxn]; clear(c, 0);
	for(RG int i = 1; i <= n; i++)
		for(RG int j = 1; j <= n; j++)
			(c[(i + j - 2) % n + 1] += 1ll * S[i] * T[j] % Mod) %= Mod;
	std::copy(c + 1, c + n + 1, S + 1);
}

int main()
{
	n = read(), s = read(), L = read(), R = read(), _x = read();
	for(RG int i = 1; i <= _x; i++) Mod *= 10;
	for(RG int i = 1; i <= n; i++) S[i] = read();
	T[1] = 1, T[2] = L, T[n] = R;
	for(; s; s >>= 1, Mul(T, T)) if(s & 1) Mul(S, T);
	for(RG int i = 1; i <= n; i++) printf("%d ", S[i]);
	return 0;
}

Extend

在這里

\(n\leq 5\times 10^4, s \leq 2^{30}\)

\(\mathrm{TLE}\)

這時我們想到，可不可以\(\mathrm{FFT}\)優化循環矩陣乘法呢？？

可以！！！

將前面的下標從\(0\)開始

\[k[x]=\sum_{(i+j) \% n = x} k[i] \times k[j] \]

這™不就是個卷積的形式嗎？？？

證明看參考博客

這里面的知識點一定要好好理(ji)解(jie lun)

這里不再贅述。

所以：

\[\begin{aligned} AB &= F\cdot diag(\hat a) \cdot F^H \cdot F \cdot diag(\hat b) \cdot F^H \\ &= F\cdot diag(\hat a) \cdot diag(\hat b) \cdot F^H \\ &= F\cdot diag(\hat a \odot \hat b) \cdot F^H \\ &= C\left( \mathcal F^{-1}(\hat a \odot \hat b)\right) \end{aligned} \]

和上面的結論是一樣的。

於是用個\(NTT\)或者\(MTT\)即可

我非常友善的\(\mathrm{mod}\;998244353\)

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#define RG register
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define clear(x, y) memset(x, y, sizeof(x))

inline int read()
{
	int data = 0, w = 1; char ch = getchar();
	while(ch != '-' && (!isdigit(ch))) ch = getchar();
	if(ch == '-') w = -1, ch = getchar();
	while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
	return data * w;
}

const int maxn(200010), Mod(998244353);
int n, s, L, R, N, r[maxn];
int S[maxn], T[maxn];

inline int fastpow(int x, int y)
{
	int ans = 1;
	while(y)
	{
		if(y & 1) ans = 1ll * ans * x % Mod;
		x = 1ll * x * x % Mod, y >>= 1;
	}
	return ans;
}

template<int opt> void FFT(int *p)
{
	for(RG int i = 0; i < N; i++) if(i < r[i]) std::swap(p[i], p[r[i]]);
	for(RG int i = 1; i < N; i <<= 1)
	{
		int rot = fastpow(3, (Mod - 1) / (i << 1));
		for(RG int j = 0; j < N; j += (i << 1))
			for(RG int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % Mod)
			{
				int x = p[j + k], y = 1ll * w * p[i + j + k] % Mod;
				p[j + k] = (x + y) % Mod, p[i + j + k] = (x - y + Mod) % Mod;
			}
	}
	if(opt == -1) std::reverse(p + 1, p + N);
}

void Mul(int *S, int *T)
{
	static int a[maxn], b[maxn];
	std::fill(a, a + N, 0); std::fill(b, b + N, 0);
	std::copy(S, S + n, a); std::copy(T, T + n, b);
	FFT<1>(a), FFT<1>(b);
	for(RG int i = 0; i < N; i++) a[i] = 1ll * a[i] * b[i] % Mod;
	FFT<-1>(a); std::fill(S, S + n, 0);
	for(RG int i = 0, inv = fastpow(N, Mod - 2); i < N; i++)
		a[i] = 1ll * a[i] * inv % Mod;
	for(RG int i = 0; i < N; i++) S[i % n] = (S[i % n] + a[i]) % Mod;
}

int main()
{
	int P = -1;
	n = read(), s = read(), L = read(), R = read();
	for(RG int i = 0; i < n; i++) S[i] = read();
	for(N = 1; N <= n + n; N <<= 1, ++P);
	for(RG int i = 0; i < N; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << P);
	T[0] = 1, T[1] = L, T[n - 1] = R;
	for(; s; s >>= 1, Mul(T, T)) if(s & 1) Mul(S, T);
	for(RG int i = 0; i < n; i++) printf("%d ", S[i]);
	return 0;
}

總結

形如\(\begin{bmatrix}c_1&c_2&\cdots&c_{n-1}&c_n\\c_2&c_3&\cdots&c_n&c_1\\\cdots&\cdots&\cdots&\cdots&\cdots\\c_n&c_1&\cdots&c_{n-2}&c_{n-1}\end{bmatrix}\)的矩陣是循環矩陣

循環矩陣的乘積仍是循環矩陣

所以我們只要維護循環矩陣的第一行

就可以\(O(n^2)\)甚至\(O(nlog_2n)\)維護循環矩陣的乘積

通過這樣一些神奇的性質，我們可以降低矩陣乘法的復雜度

來出毒瘤題

沒了