Cayley-Hamilton 定理学习笔记
Cayley-Hamilton 定理学习笔记
数学好难
Cayley-Hamilton 定理是一个家喻户晓的定理，可以用来加速矩阵快速幂。一般来说可以加速到 \(O(k^4+k^2\log n)\) 的时间复杂度，据说可以加速到 \(O(k^3+k\log k\log n)\)，但是我不会（
Cayley-Hamilton 定理
Cayley-Hamilton 定理即：设 \(A\) 是一个 \(k\) 阶矩阵，则其特征多项式为 \(f(x)=|xE-A|=x^k+c_1x^{k-1}+c_2x^{k-2}+\cdots+c_{k-1}x+c_k\)，其中 \(E\) 是单位矩阵，\(c\) 是一系列系数，\(x\) 可以取很多种值，如其可以取为一个矩阵或一个数，注意当它取成一个矩阵的时候 \(f(x)\) 也是一个矩阵（尽管 \(|xE-A|\) 看上去是一个数，具体为什么是这样的小编也不清楚）。
我们发现，\(f(A)=|AE-A|=0\)，也即 \(A^k+c_1A^{k-1}+c_2A^{k-2}+\cdots+c_{k-1}A+c_k=0\)。
求矩阵特征多项式
那我咋求一个矩阵的特征多项式呢？
从特征多项式的名字可以看出它是个多项式（），因此我们可以考虑插值。不难发现它是个 \(k\) 次多项式，因此我们只要带 \(k+1\) 个点值就行了。
再看看上面的式子，我们发现特征多项式其实就是一个行列式。因此我们可以随便带 \(k+1\) 个值，对每个值暴力高斯消元，时间复杂度 \(O(k^4)\)。
有通过构造相似上海森堡矩阵使这一步加速到 \(O(k^3)\) 的方法，但是我不会。这篇博客里面有讲。
加速矩阵快速幂
假设我们要求 \(M^n\)，其中 \(M\) 是一矩阵。
首先肯定得先求出 \(f()\)，这里时间是 \(O(k^4)\)。
我们知道 \(f(M)=|ME-M|=0\)，因此有 \(M^n=(M^n\bmod f(M))\)。此处取模类似于多项式取模。
换种表述方法，设 \(f(x)\) 是 \(M\) 的特征多项式（与上文定义相同），\(g(x)=x^n\)，则有 \(BM^n=B\cdot g(M)\)。设 \(g(x)=f(x)h(x)+R(x)\)，因为 \(f(M)=0\)，所以 \(g(M) = R(M)\)，也即 \(g(m)=(g\bmod f)(M)\)。
考虑平时做模意义下整数快速幂的方法，将整数乘法换成多项式乘法，整数取模换成多项式取模，即可求出 \((g\bmod f)(x)\)。如果暴力乘和取模这一步是 \(k^2\log n\) 的，换成 NTT/FFT 和常用的多项式取模即可做到 \(O(k\log k\log n)\)。
最后得把每个 \(M^i\) 给带进去，一共要带 \(k\) 个，因为 \(f()\) 是 \(k\) 次多项式。如果能通过一些奇技淫巧把每个 \(M^i\) 搞出来的话时间复杂度是 \(O(k^3)\)，暴力矩乘则是 \(O(k^4)\)。
最终复杂度 \(O(k^4+k^2\log n)\)。这并不是理论最优复杂度，但是再快的我就不会了（
加速齐次线性递推
所有齐次线性递推的转移矩阵的形态都很相似，因此其特征多项式形态也都很相似。
假设一数列 \(h\) 满足 \(h_n=\sum_{i=1}^ka_ih_{n-i}\)，\(n\ge k\)。
我们假设转移矩阵为 \(M\)，初始的列向量为 \(B=\begin{bmatrix}h_0\\h_1\\ \vdots\\ h_{k-1}\end{bmatrix}\)，并设我们要求 \(h_n\)，那么我们要求的即是 \(BM^n[0]\) 。
同上，设 \(f(x)\) 是 \(M\) 的特征多项式（与上文定义相同），\(g(x)=x^n\)，则 \(BM^n=B\cdot(g\bmod f)(M)\)。
那我们咋求 \((g\bmod f)(M)\) 呢？
首先我们肯定得先把 \(f(x)\) 求出来。
不妨将转移矩阵写出来，我们要求的即是
也不是不能直接暴力 \(O(k^4)\) 求，但是一是当 \(k\) 略大的时候就挂了，二是直接暴力求有点对不起形态这么优美的矩阵（
不妨设 \(F_u(x)\) 为 取其最后 \(u\) 行 \(u\) 列的主子式的特征多项式，不妨设其为 \(\sum_{i=0}^{u}c_{u,i}x^i\)。
由行列式定义式，矩阵第一行取 \(x\) 的贡献与取 \(-1\) 的贡献加起来就是整个行列式。
矩阵第一行取 \(x\) 时，我们发现此时第一行与第一列被除去了，其余没有变化，于是贡献为 \(xF_{k-1}(x)\)。
矩阵第一行取 \(-1\) 时，我们发现此时第二行也只能选 \(-1\)，然后第三行也只能选 \(-1\)，以此类推，前 \(k-1\) 行都只能选 \(-1\)，于是最后一行只能选 \(-a_k\)。在这里面，第一行和最后一行是一对逆序对，第二行和最后一行是一对逆序对，……，共有 \(n-1\) 对逆序对，于是贡献为 \((-1)^ka_k\cdot(-1)^{k-1}=-a_k\)。
由此，我们发现：
好家伙，原来知道了 \(a\) 就知道了 \(f(x)\)。
再回过头来看看我们要做什么：我们要求 \(BM^n[0]\)。\(M^n\) 可以直接用上文方法求，求出来的式子大概可以表示为 \(\sum_{i=0}^{k-1}c_iM^i\)，于是 \(BM^n=B\sum_{i=0}^{k-1}c_iM^i=\sum_{i=0}^{k-1}c_iBM^i\)。
我们发现，\(BM^i[0]=B[i]\)，又上式中有 \(i < k\)，因此只要知道了 \(h_0,h_1,\cdots,h_{k-1}\) 也就知道了 \(BM^n[0]\)。同理，要是想知道 \(BM^n\) 整个列向量，我们也只需知道 \(h_0,h_1,\cdots,h_{2k-1}\)。
若只求单个元素，总时间复杂度是求 \(M^n\) 时的 \(O(k^2\log n)/O(k\log k\log n)\)。若要求整个列向量，要么直接用 \(\sum_{i=0}^{k-1}c_iBM^i\) 这个式子暴力矩阵乘法，时间复杂度 \(O(k^2\log n+k^3)\)，比直接矩阵快速幂少个 \(\log\) 或者 \(k\)，且这样我们就只需要知道 \(h_0,h_1,\cdots,h_{k-1}\)；要么挨个元素去算，时间复杂度 \(O(k^3\log n)/O(k^2\log k\log n)\)。可以发现不用 NTT 和取模的时间和矩快一样，用了的则是把 \(k\) 变成了 \(\log k\)。
还是搞个例题免得字数太少：BZOJ4161
板题。\(h\) 和 \(a\) 都已经知道了，直接套上述方法就行了。
自己写的傻逼代码：
#include <bits/stdc++.h>
using namespace std;
int n , k;
long long ans;
struct vec
{
long long a[4400];
void clear()
{
memset(a , 0 , sizeof(a));
}
} a , b , e , P;
const int mod = 1e9 + 7;
vec mul( vec a , vec b , vec P )
{
vec ans; ans.clear();
for(int i = 0 ; i <= k ; i++ )
for(int j = 0 ; j <= k ; j++ )
(ans.a[i + j] += a.a[i] * b.a[j] % mod) %= mod;
for(int i = 2 * k ; i >= k ; i-- )
{
if(ans.a[i])
{
int qwq = ans.a[i];
for(int j = 0 ; j <= k ; j++ )
ans.a[i - j] = (ans.a[i - j] - qwq * P.a[k - j] % mod + mod) % mod;
}
}
return ans;
}
vec exp( vec a , int b , vec P )
{
vec ans = e , t = a;
while(b)
{
if(b & 1) ans = mul(ans , t , P);
t = mul(t , t , P); b >>= 1;
}
return ans;
}
int main()
{
e.a[0] = 1;
scanf("%d%d" , &n , &k);
for(int i = 1 ; i <= k ; i++ )
{
scanf("%lld" , &a.a[i]);
P.a[k - i] = a.a[i] = (mod - a.a[i]) % mod;
}
for(int i = 0 ; i < k ; i++ ) scanf("%lld" , &b.a[i]) , b.a[i] = (b.a[i] + mod) % mod;
P.a[k] = 1;
a.clear(); a.a[1] = 1;
a = exp(a , n , P);
for(int i = 0 ; i < k ; i++ ) ans = (ans + b.a[i] * a.a[i] % mod) % mod;
printf("%lld" , ans);
return 0;
}
/*
*/
总算把这玩意学懂了...可能过不了两天又要忘了XDD
