矩阵快速幂——斐波那契数列进一步优化

快速幂优化矩阵幂、乘法

---

对于一般的矩阵计算有 \(A_{m,n}*B_{n,p}=C_{m,p}\)，其中作为乘积因子的两个矩阵必须满足前因子列数与后因子行数相同

积的行数等于前因子的行数，列数等于后因子的列数 ，任意的 \(c_{i,j}\) 可由定义的计算得出 \(c_{i,j}=\sum_{k=1}^{n} a_{i,k}*b_{k,j}\)

---

幂运算

可以通过快速幂优化矩阵的幂乘法，即计算 \(A^k\)

struct matrix {
    long long mtx[N][N];
    matrix() { memset(mtx, 0, sizeof mtx); }
};//定义矩阵结构体

其中 matrix() { memset(mtx, 0, sizeof mtx); }可以在每次定义矩阵时，都对矩阵成员进行初始化为 \(0\)

matrix operator*(matrix& x, matrix& y) {
    matrix t;
    for (int i = 1; i <= n; i++)//行
        for (int j = 1; j <= n; j++)//列
            for (int k = 1; k <= n; k++)//当前矩阵的内数
                t.mtx[i][j] += x.mtx[i][k] * y.mtx[k][j];
    return t;
}

重载矩阵的乘法运算符，通过普遍公式计算给出

matrix matrix_quick_pow(matrix a, long long k) {//矩阵类型的快速幂计算函数
    matrix res;
    for (int i = 1; i <= n; i++) res.mtx[i][i] = 1;
    while (k) {
        if (k & 1) res = res * a;
        a = a * a;
        k >>= 1;
    }
    return res;
}

定义 res 时，需要对其进行初始化为一个单位矩阵作为乘积的基

---

加速斐波那契数列

对于斐波那契数列有以下通项公式 \(F_n=F_{n-1}+F_{n-2}\)

设存在矩阵 \(A\) 使得

\[\left[\begin{matrix}F_n&F_{n-1}\end{matrix}\right] = \left[\begin{matrix}F_{n-1}&F_{n-2}\end{matrix}\right] * A \quad A = \left[\begin{matrix}a&b\\c&d\end{matrix}\right] \]

结合通项公式，合并矩阵后等式变化为

\[\left[\begin{matrix}F_n&F_{n-1}\end{matrix}\right] = \left[\begin{matrix}F_{n-1}&F_{n-2}\end{matrix}\right] * \left[\begin{matrix}1&1\\1&0\end{matrix}\right] \]

因为\(F_1=F_2=1\) 已知，递推到底层有

\[\left[\begin{matrix}F_n&F_{n-1}\end{matrix}\right] = \left[\begin{matrix}F_{2}&F_{1}\end{matrix}\right] * \left[\begin{matrix}1&1\\1&0\end{matrix}\right]^{n-2} \]

最终可化简为

\[\begin{cases} \left[\begin{matrix}F_n&F_{n-1}\\F_{n-1}&F_{n-2}\end{matrix}\right] = \left[\begin{matrix}1&1\\1&0\end{matrix}\right]^{n-1}&(n>2) \\\\F_n=1&(n=1,2) \end{cases} \]

此时可以发现，计算斐波那契数列的各项值就只需要计算矩阵的 \(n\) 次方，结合上面的矩阵快速幂，可以做到在 \(O(\log n)\) 的时间复杂度内计算出第 \(n\) 项

---

matrix a;
init(a);
void init(matrix &a) {
	a.mxt[1][1] = 1;
	a.mxt[1][2] = 1;
	a.mxt[2][1] = 1;
}

初始化斐波那契矩阵

auto ans = matrix_quick_pow(a, n - 1);
cout << ans.mxt[1][1] << endl;

根据公式，最后输出第一行第一列的元素即可

完整代码：

#include <iostream>
#include <cstring>
using namespace std;

const long long mod = 1e9 + 7;
long long n;

struct matrix {
	long long mxt[3][3];
	matrix() { memset(mxt, 0, sizeof mxt); };
};

matrix operator*(matrix& x, matrix& y) {
	matrix t;
	for (int i = 1; i <= 2; i++)
		for (int j = 1; j <= 2; j++)
			for (int k = 1; k <= 2; k++)
				t.mxt[i][j] = (t.mxt[i][j] + x.mxt[i][k] * y.mxt[k][j]) % mod;
	return t;
}

void init(matrix &a) {
	a.mxt[1][1] = 1;
	a.mxt[1][2] = 1;
	a.mxt[2][1] = 1;
}

matrix matrix_quick_pow(matrix a, long long n) {
	matrix res;
	for (int i = 0; i <= 2; i++)
		res.mxt[i][i] = 1;
	while (n)
	{
		if (n & 1) res = res * a;
		a = a * a;
		n >>= 1;
	}
	return res;
}

int main()
{
	matrix a;
	init(a);
	cin >> n;
	auto ans = matrix_quick_pow(a, n - 1);
	cout << ans.mxt[1][1] << endl;
	return 0;
}