矩阵快速幂的常数优化——对角化与 Jordan 分解
只会无脑写出 \(5 \times 5\)（或者更大）的转移矩阵而不会推式子？还在尝试卡常？本文提供了一种方法在这些情况下优化常数从而能够强行 AC。
（其实等价于推式子，但更无脑。）
对角化
理论
给定 \(n \times n\) 方阵 \(A\)，若其可对角化，则可分解为 \(A = PDP^{-1}\)，其中 \(D\) 为对角矩阵，\(P\) 可逆。于是 \(A^{k} = P D^k P^{-1}\)，只需计算 \(D^k\)，而对角矩阵的幂无需 \(O(n^3 \log k)\) 的矩阵快速幂，只要将 \(O(n)\) 个对角元素分别快速幂即可。于是我们就将 \(O(n^3 \log k)\) 优化到了 \(O(n \log k)\)。
实践
假设你想求出 Fibonacci 数列第 \(n\) 项 \(Fib_n\) 的值。写出矩阵
\[A = \begin{bmatrix}
1 & 1 \\
1 & 0
\end{bmatrix}\]
后可以用 Python 的 SymPy 库来计算。
import sympy as sp
A = sp.Matrix([
[1, 1],
[1, 0]
])
P, D = A.diagonalize()
P.simplify()
D.simplify()
P_inv = P.inv()
P_inv.simplify()
sp.pprint(P)
sp.pprint(D)
sp.pprint(P_inv)
解得
\[P=\begin{bmatrix}
\dfrac{1-\sqrt5}{2} & \dfrac{1+\sqrt5}{2} \\
\\ 1 & 1
\end{bmatrix}
\\
D=\begin{bmatrix}
\dfrac{1-\sqrt5}{2} & 0 \\
0 & \dfrac{1+\sqrt5}{2}
\end{bmatrix}\\
P^{-1} = \begin{bmatrix}
-\dfrac{\sqrt5}{5} & \dfrac{\sqrt5}{10}+\dfrac{1}{2} \\
\\ \dfrac{\sqrt5}{5} & \dfrac{1}{2} - \dfrac{\sqrt5}{10}
\end{bmatrix}
\]
于是
\[\begin{aligned}
& A^n \\
=& P D^n P^{-1} \\
=& \begin{bmatrix}
\dfrac{1-\sqrt5}{2} & \dfrac{1+\sqrt5}{2} \\
\\ 1 & 1
\end{bmatrix}
\begin{bmatrix}
\left(\dfrac{1-\sqrt5}{2}\right)^n & 0 \\
0 & \left(\dfrac{1+\sqrt5}{2}\right)^n
\end{bmatrix}
\begin{bmatrix}
-\dfrac{\sqrt5}{5} & \dfrac{\sqrt5}{10}+\dfrac{1}{2} \\
\\ \dfrac{\sqrt5}{5} & \dfrac{1}{2} - \dfrac{\sqrt5}{10}
\end{bmatrix}
\end{aligned}
\]
于是只要求两个元素的快速幂即可。在取模时若 \(5\) 不为模意义下二次剩余，可以直接无脑扩域。
事实上，将这个式子再化简一下就可以得到 Fibonacci 数列的通项公式：
\[Fib_n = \dfrac{1}{\sqrt5}\left( \left( \dfrac{1+\sqrt5}{2} \right)^n - \left( \dfrac{1-\sqrt5}{2} \right)^n \right)
\]
Jordan 分解
理论
如果 \(A\) 不能对角化呢？类似地，我们有 Jordan 分解，可以将矩阵分解为 \(A=PJP^{-1}\)，其中 \(J\) 是 Jordan 矩阵，\(P\) 可逆。
Jordan 矩阵是由 Jordan 块组成的对角矩阵。Jordan 块可以写成对角矩阵与幂零矩阵的和，因此可以用二项式定理展开并截断到幂零矩阵不为零的部分。
实践
以 lg4834 为例。题目要求
\[\dfrac{2}{n(n+1)} \sum \limits_{i=1}^{n}i \mathrm{Fib}(i)
\]
对 \(998244353\) 取模后的值。可以无脑设状态为
\[\mathbf{v} = \begin{bmatrix}
\mathrm{Fib}(n), \mathrm{Fib}(n-1), n\mathrm{Fib}(n), (n-1) \mathrm{Fib}(n-1), \sum \limits_{i=1}^{n-1}i\mathrm{Fib}(i)
\end{bmatrix}^T
\]
并写出转移矩阵
\[A=\begin{bmatrix}
1 & 1 & 0 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 \\
1 & 2 & 1 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 & 1
\end{bmatrix}
\]
但是矩阵太大，直接做会 TLE，并且它也无法对角化。这时可以用 SymPy 求它的 Jordan 分解：
import sympy as sp
A = sp.Matrix([
[1, 1, 0, 0, 0],
[1, 0, 0, 0, 0],
[1, 2, 1, 1, 0],
[0, 0, 1, 0, 0],
[0, 0, 1, 0, 1]
])
P, J = a.jordan_form()
P.simplify()
J.simplify()
P_inv = P.inv()
P_inv.simplify()
sp.pprint(P)
sp.pprint(J)
sp.pprint(P_inv)
发现
\[J = \begin{bmatrix}
1& 0& 0& 0& 0 \\
0&\dfrac{1-\sqrt5}{2}& 1& 0& 0 \\
0& 0&\dfrac{1-\sqrt5}{2}& 0& 0 \\
0& 0& 0& \dfrac{1+\sqrt5}{2}& 1 \\
0& 0& 0& 0& \dfrac{1+\sqrt5}{2}
\end{bmatrix}
\]
其对应的幂零矩阵
\[\begin{bmatrix}
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 0
\end{bmatrix}
\]
的幂零指数为 \(2\)，因此
\[J^k = \begin{bmatrix}
1 & 0 & 0 & 0 & 0 \\
0 & \lambda_1^k & k \lambda_1^{k-1} & 0 & 0 \\
0 & 0 & \lambda_1^k & 0 & 0 \\
0 & 0 & 0 & \lambda_2^k & k \lambda_2^{k-1} \\
0 & 0 & 0 & 0 & \lambda_2^k
\end{bmatrix}
\]
其中 \(\lambda_1 = \dfrac{1-\sqrt5}{2}, \lambda_2 = \dfrac{1+\sqrt5}{2}\)。
本来矩阵乘法有 \(125\) 的常数，现在优化后只要做 \(2\) 次快速幂即可。
结果里出现了 \(\sqrt5\)，在 \(\bmod \ 998244353\) 意义下需要扩域。
struct CC{ // Z_998244353[sqrt(5)]
int re, im; // i = sqrt(5)
CC(){}
constexpr CC(int x, int y=0): re(x), im(y){}
CC &operator+=(CC b){ re += b.re; im += b.im; re %= MOD; im %= MOD; return *this; }
CC operator+(CC b) const { return b += *this; }
CC operator-() const { return CC(MOD-re, MOD-im); }
CC &operator-=(CC b){ return *this += -b; }
CC operator-(CC b) const { return -b += *this; }
CC operator*(CC b) const {
ll ree = 1ll * re * b.re + 5ll * im * b.im;
ll imm = 1ll * re * b.im + 1ll * im * b.re;
return CC(ree % MOD, imm % MOD);
}
CC &operator*=(CC b) { return *this = *this * b; }
};
矩阵运算与定义部分：
struct Matrix55{
CC elem[5][5];
Matrix55(){}
Matrix55(CC val[]){
UP(i, 0, 5) UP(j, 0, 5) elem[i][j] = val[i*5+j];
}
Matrix55 operator*(Matrix55 const& b){
Matrix55 tmp;
UP(i, 0, 5) UP(j, 0, 5) tmp.elem[i][j] = 0;
UP(i, 0, 5) UP(j, 0, 5) UP(k, 0, 5){
tmp.elem[i][k] += elem[i][j] * b.elem[j][k];
}
return tmp;
}
};
constexpr CC inv2 = 499122177, inv5 = 598946612, inv10 = 299473306, sqrt5 = CC(0, 1);
Matrix55 Mat_P = [](){
// Matrix([
// [0, 0, 1, 0, 1],
// [0, 0, -sqrt(5)/2 - 1/2, 0, -1/2 + sqrt(5)/2],
// [0, 1/2 - sqrt(5)/2, 3/2 - sqrt(5)/2, 1/2 + sqrt(5)/2, sqrt(5)/2 + 3/2],
// [0, 1, 1, 1, 1],
// [1, 3/2 - sqrt(5)/2, 0, sqrt(5)/2 + 3/2, 0]])
CC val[] = {
0, 0, 1, 0, 1,
0, 0, -inv2 - sqrt5 * inv2, 0, -inv2 + sqrt5 * inv2,
0, inv2 - sqrt5 * inv2, inv2 * 3 - sqrt5 * inv2, inv2 + sqrt5 * inv2, inv2 * 3 + sqrt5 * inv2,
0, 1, 1, 1, 1,
1, inv2 * 3 - sqrt5 * inv2, 0, sqrt5 * inv2 + inv2 * 3, 0
};
return Matrix55(val);
}();
Matrix55 Mat_P_inv = [](){
// Matrix([
// [ 3, 1, -1, -1, 1],
// [-1/2 + 3*sqrt(5)/10, sqrt(5)/5, -sqrt(5)/5, sqrt(5)/10 + 1/2, 0],
// [ 1/2 - sqrt(5)/10, -sqrt(5)/5, 0, 0, 0],
// [-3*sqrt(5)/10 - 1/2, -sqrt(5)/5, sqrt(5)/5, 1/2 - sqrt(5)/10, 0],
// [ sqrt(5)/10 + 1/2, sqrt(5)/5, 0, 0, 0]])
CC val[] = {
3, 1, MOD-1, MOD-1, 1,
-inv2 + sqrt5*inv10*3, sqrt5*inv5, -sqrt5*inv5, sqrt5*inv10 + inv2, 0,
inv2 - sqrt5*inv10, -sqrt5*inv5, 0, 0, 0,
-sqrt5*inv10*3 - inv2, -sqrt5*inv5, sqrt5*inv5, inv2 - sqrt5*inv10, 0,
sqrt5*inv10 + inv2, sqrt5*inv5, 0, 0, 0
};
return Matrix55(val);
}();
std::function<Matrix55(int)> calc_J_n = [](int n){
CC lam1 = inv2 - sqrt5 * inv2, lam2 = inv2 + sqrt5 * inv2;
CC lam1_n = qpow(lam1, n-1), lam2_n = qpow(lam2, n-1);
CC val[] = {
1, 0, 0, 0, 0,
0, lam1_n * lam1, lam1_n * n, 0, 0,
0, 0, lam1_n * lam1, 0, 0,
0, 0, 0, lam2_n * lam2, lam2_n * n,
0, 0, 0, 0, lam2_n * lam2
};
return Matrix55(val);
};
单次询问时注意特判 \(n \le 2\) 的情况：
int in;
void work(){
cin >> in;
if(in <= 2){ cout << 1 << '\n'; return; }
CC x = qpow(CC(in+1) * in, MOD-2) * 2;
auto p = Mat_P * calc_J_n(in-1) * Mat_P_inv;
CC s = p.elem[4][0] + p.elem[4][1] + p.elem[4][2]*2 + p.elem[4][3] + p.elem[4][4];
s *= x;
cout << s.re << '\n';
}
于是我们在不推结论的情况下，用 \(5 \times 5\) 的矩阵乘法通过了此题。评测记录
浙公网安备 33010602011771号