矩阵快速幂的常数优化——对角化与 Jordan 分解
cnblogs
本文同时作为洛谷 P4834 的题解（然而洛谷只能投一种分类）。
只会无脑写出
（其实等价于推式子，但更无脑。）
对角化
理论
给定
实践
假设你想求出 Fibonacci 数列第
后可以用 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)解得
于是
于是只要求两个元素的快速幂即可。在取模时若
事实上，将这个式子再化简一下就可以得到 Fibonacci 数列的通项公式：
Jordan 分解
理论
如果
Jordan 矩阵是由 Jordan 块组成的对角矩阵。Jordan 块可以写成对角矩阵与幂零矩阵的和，因此可以用二项式定理展开并截断到幂零矩阵不为零的部分。
实践
以 lg4834 为例。题目要求
对
并写出转移矩阵
但是矩阵太大，直接做会 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)发现
其对应的幂零矩阵
的幂零指数为
其中
本来矩阵乘法有
结果里出现了
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);
};单次询问时注意特判
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';
}于是我们在不推结论的情况下，用