矩陣

定義

在數學中，矩陣（Matrix）是一個按照長方陣列排列的復數或實數集合。

就像下面一樣：

\[\begin{bmatrix} 1&2&3\\1&2&3 \end{bmatrix}\]

特別的，主對角線上為 1，其余位置為 0 的叫做單位矩陣 \(I\)：

\[\begin{bmatrix} 1&0&\cdots&0\\ \vdots&1&\ddots&0\\ 0&0&\dots&1\\ \end{bmatrix}\]

矩陣的運算

矩陣的加減法是逐個元素進行的。

乘法

矩陣相乘只有在第一個矩陣的列數和第二個矩陣的行數相同時才有意義。

設 \(P\) 為 \(M \times N\) 的矩陣，\(Q\) 是 \(N \times Q\) 的矩陣，\(C\) 矩陣為 \(C = PQ\)，則 \(C\) 矩陣上 \(i\) 行 \(j\) 列的元素為。

\[c_{i,j} = \sum\limits_{k=1}^{M}P_{i,k}Q_{k,j} \]

沒有看懂，沒事，有一種比較簡單的方法。

設 \(A\) 矩陣為:

\[\begin{bmatrix} a&b\\ c&d\\ \end{bmatrix}\]

設 \(B\) 矩陣為:

\[\begin{bmatrix} e\\f \end{bmatrix}\]

所以 \(A\) 與 \(B\) 的乘積是：

\[e\begin{bmatrix} a\\c \end{bmatrix} + f\begin{bmatrix} b\\d \end{bmatrix}\]

注意，矩陣乘法滿足結合律，不滿足一般的交換律。

給出一個封裝好的矩陣乘法：

struct mat{
	long long c[N][N];
    void unit(){
        memset(c,0,sizeof(c));
        for(int i=1;i<=n;i++)   c[i][i] = 1;
    }
    void mat(){memset(c,0,sizeof(c));}
	inline mat operator*(const mat &M){
		mat ans;
		memset(ans.c,0,sizeof(ans.c);
		for(int i=1;i<=n;i++)
			for(int j=1;j<=n;j++)
				for(int k=1;k<=n;k++)
					ans.c[i][j] += (c[i][k] * M.c[k][j]) % P,ans.c[i][j] %= P;
		return ans;
	}
};

如果算一個矩陣的 \(A^k\)，可以使用快速冪來加速到 \(\mathcal{O}(n^3 \log k)\)。

// 快速冪
mat qpow(mat x,int y){
    mat ans.unit();
    while(y > 0){
        if(y % 2 == 1)
            ans = ans * x;
        y = y / 2;
        x = x * x;
    }
    return ans;
}

矩陣的應用

快速求斐波拉契數列

\(F_1 = F_2 = 1\)，\(F_i = F_{i-1} + F_{i-2} \ (i \geq 3)\) 就是斐波拉契數列。

當然可以使用遞推的方法來求，時間復雜度 \(\mathcal{O}(n)\) ，還挺快。

但是遇到這道題目就不行了，\(2^{63}\) 絕對會 TLE，所以可以用矩陣加速遞推。

我們可以設 \(f(n) = \begin{bmatrix}F_{n} & F_{n-1}\end{bmatrix}\) ，我們希望可以從 \(f(n-1)\) 推出它。

正好我們之前學了矩陣的乘法，所以我們可以嘗試設一個矩陣 \(\text{base}\) ，使得 \(f(n) = \text{base} \times f(n-1)\)，也就是：

\[\begin{bmatrix} F_i & F_{i-1} \end{bmatrix} = \text{base} \times \begin{bmatrix} F_{i-1} & F_{i-2} \end{bmatrix}\]

因為前后都是 \(1 \times 2\) 的矩陣，所以可以知道 \(\text{base}\) 矩陣是 \(2 \times 2\) 的矩陣。

所以 \(\text{base} = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix}\)

所以等式轉化為：

\[\begin{bmatrix} F_i & F_{i-1} \end{bmatrix} = F_{i-1} \times \begin{bmatrix} a & b \end{bmatrix} + F_{i-2} \times \begin{bmatrix} c &d \end{bmatrix}\]

\[\begin{bmatrix} F_i & F_{i-1} \end{bmatrix} = \begin{bmatrix} a F_{i-1} + c F_{i-2} & b F_{i-1} + d F_{i-2} \end{bmatrix}\]

又因為 \(F_i = F_{i-1} + F_{i-2}\)，所以 \(a = 1,c = 1,b = 1,d = 0\)。

綜上所述， \(\text{base} = \begin{bmatrix} 1 & 1\\ 1 & 0\end{bmatrix}\)。

那么因為初始項為 \(F_1,F_2\)。所以 \(F_n = \begin{bmatrix} F_1 & F_2 \end{bmatrix} \times \text{base}^{n-2}\)

好好品品。

\(\text{base}^{n-2}\) 矩陣可以通過上面的快速冪來求，所以時間復雜度為 \(\mathcal{O}(\log n)\)。

這題核心代碼：

struct mat{
	long long c[4][4];
	mat() {memset(c,0,sizeof(c));}
	inline mat operator*(const mat& M){
		mat res;
		for(int i=1;i<=2;i++)
			for(int j=1;j<=2;j++)
				for(int k=1;k<=2;k++)
					res.c[i][j] = (res.c[i][j] + c[i][k]*M.c[k][j]) % mod;
		return res;
	}
};

mat base,ans;

void init(){
	base.c[1][1] = base.c[1][2] = base.c[2][1] = 1;
	ans.c[1][1] = ans.c[1][2] = 1; 
}

void qpow(long long y){
	while(y > 0){
		if(y % 2 == 1)
			ans = ans * base;
		y = y / 2;
		base = base * base;
	}
}

其實矩陣可以干事情的很多，就寫這么多吧