浅谈循环矩阵的乘法和快速幂

前言

前置技能：矩阵乘法，矩阵快速幂

当然你不会的话也不会点进来（滑稽）

今天上午的$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$都是循环矩阵，那么：

$a+b$是一个循环矩阵

$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;
}