1. 高斯消元
定义以下三种操作为 初等行变换 ：
- 将某一行乘上 (c)
- 将某一行的 (c) 倍加到另一行
- 交换两行
这里介绍高斯-约旦消元法，利用初等行变换将矩阵转化成对角矩阵，然后就很容易求出解。
- 枚举 (i) 表示现在要消去第 (i) 项
- 在 ([i,n]) 内选择一个 (k) 满足 (a[k][i] eq 0)，若不存在则说明无解
- 交换第 (i) 行和第 (k) 行
- 对于 ([i+1,n]) 内所有 (j)，将第 (j) 行减去 (dfrac{a[j][i]}{a[i][i]}) 倍的第 (i) 行
容易发现，最后第 (i) 项系数只会出现在第 (i) 行。并且任意时刻，第 (i) 行的 (1dots i-1) 项系数都为 (0)，所以不会对其他方程的这些项造成影响。
inline bool zero(double t) {
return fabs(t) < 1e-8;
}
signed main() {
scanf("%d", &n);
rep(i, 1, n) rep(j, 1, n + 1)
scanf("%lf", &a[i][j]);
rep(i, 1, n) {
k = i;
if(zero(a[i][i])) rep(j, i + 1, n) if(!zero(a[j][i])) { k = j; break; }
if(i ^ k) rep(j, i, n + 1) swap(a[i][j], a[k][j]);
if(zero(a[i][i])) return puts("No Solution") & 0;
rep(j, 1, n) if(i != j) {
t = a[j][i] / a[i][i];
rep(k, i, n + 1) a[j][k] -= a[i][k] * t;
}
}
rep(i, 1, n) printf("%.2lf
", a[i][n + 1] / a[i][i]);
return 0;
}
2. 行列式
对于一个方阵，其行列式是线性变换的伸缩因子。例如在三维空间中，行列式表示线性变换造成的体积放大倍率。
代数定义：(det A=|A|=sumlimits_{p是n的排列}(-1)^{p的逆序对数}prodlimits_{i=1}^na_{i,p_i})
看右边那串东西，其实是从方阵中选出 (n) 个数，使得没有两个数在同一行或同一列，然后把它们乘起来。所以如果选出的数中有一个 (0) 就没有贡献了。
排列的性质
- 交换两数，逆序对数的奇偶性改变
- (nge 2) 时，奇偶排列（逆序对数为奇 / 偶数的排列）各占一半
行列式的性质
- (|M|=|M^T|)
- 交换两行，行列式反号
- 用一个数乘行列式的某行，等于用这个数乘此行列式
- 如果行列式中有两行成比例，则行列式等于 (0)（先把比例系数提出来，发现两行相同了，交换两行行列式不变。但根据性质 2，行列式反号了，所以只能为 (0)）
- 若行列式中某一行是两组数之和，则这个行列式等于两个行列式之和，这两个行列式分别以这两组数为该行，而其余各行与原行列式对应各行相同。即：(left|egin{matrix}dots\a_{i,1}+b_{i,1},a_{i,2}+b_{i,2}dots a_{i,n}+b_{i,n}\dotsend{matrix} ight|=left|egin{matrix}dots\a_{i,1},a_{i,2}dots a_{i,n}\dotsend{matrix} ight|+left|egin{matrix}dots\b_{i,1},b_{i,2}dots b_{i,n}\dotsend{matrix} ight|)
- 将某一行的倍数加到另一行，行列式的值不变
- (left|egin{matrix}A~0\C~Bend{matrix} ight|=|A||B|)，其中 (A,B) 是方阵 （从选数的角度考虑）
- 一个矩阵的行列式可以按行展开，即：(left|egin{matrix}dots\a_{i,1},a_{i,2}dots a_{i,n}\dotsend{matrix} ight|=left|egin{matrix}dots\a_{i,1},0dots 0\dotsend{matrix} ight|+left|egin{matrix}dots\0,a_{i,2},0dots 0\dotsend{matrix} ight|+dots+left|egin{matrix}dots\0,dots 0,a_{i,n}\dotsend{matrix} ight|)（还是从选数的角度理解，第 (i) 行只能选一个数，选 (a_{i,1},a_{i,2},dots,a_{i,n}) 的方案是独立的）。进一步地，定义代数余子式：(A_{i,j}=(-1)^{i+j}|M_{i,j}|)，其中 (M_{i,j}) 是原矩阵删去第 (i) 行和第 (j) 列得到的矩阵。一个矩阵的行列式，等于它某一行中，每个元素乘上它的代数余子式的和。即：(|B|=sumlimits_{j=1}^nb_{i,j}A_{i,j}=(-1)^{i+j}b_{i,j}|M_{i,j}|)（展开后，每个矩阵的第 (i) 行只有一个数，用 (i+j-2) 次交换邻行邻列把它交换到 ((1,1))，然后用性质 7）
- 上（下）三角矩阵的行列式为其主对角线元素之积。
行列式求值
现在考虑如何计算一个矩阵的行列式。观察到第 9 条性质，我们可以把矩阵消成上三角矩阵，同时注意符号的变化，然后直接得出答案。
但当矩阵内都是整数时，仍会涉及实数运算，导致精度较低。考虑一种辗转相除的方法：每次选定两行 (i,j)，将第 (j) 行尽可能多地减去第 (i) 行，然后交换这两行，直到 (a_{i,i}=0)。这样时间复杂度是 (O(n^3+n^2log V)) 的。
rep(i, 1, n) rep(j, i + 1, n) {
while(a[i][i]) {
t = a[j][i] / a[i][i];
rep(k, i, n) a[j][k] = (a[j][k] - a[i][k] * t % p + p) % p;
swap(a[i], a[j]);
ans = -ans;
}
swap(a[i], a[j]);
ans = -ans;
}
if(ans < 0) ans += p;
rep(i, 1, n) (ans *= a[i][i]) %= p;
3. 矩阵求逆
对于矩阵 (A)，若存在矩阵 (B) 使得 (AB=I)（(I) 是单位矩阵），则 (B) 是 (A) 的逆，也可写作 (A^{-1})。
对矩阵进行高斯消元时，我们发现可以仅用初等行变换将矩阵变成单位矩阵，而每一种初等行变换都可以表示成一次矩阵乘法。即 (AP_1P_2dots P_k=I)。那么 (P_1P_2dots P_k) 不就是答案吗？
我们可以在高斯消元时记录下所有操作，然后搞一个单位矩阵，再把这些操作做一遍。或者更方便地，我们在原矩阵右边添加一个单位矩阵，把左边消成单位矩阵，右边就求出了逆。
若原矩阵不能消成单位矩阵，则原矩阵不可逆。
n = read();
rep(i, 1, n) rep(j, 1, n) a[i][j] = read();
rep(i, 1, n) a[i][i + n] = 1;
rep(i, 1, n) {
if(!a[i][i]) rep(j, i + 1, n) if(a[j][i]) { swap(a[i], a[j]); break; }
if(!a[i][i]) return puts("No Solution") & 0;
tmp = poww(a[i][i]);
rep(j, 1, n) if(i != j) {
t = a[j][i] * tmp % P;
rep(k, 1, n * 2) a[j][k] = (a[j][k] - a[i][k] * t % P + P) % P;
}
rep(j, i, n * 2) (a[i][j] *= tmp) %= P;
}
rep(i, 1, n) {
rep(j, n + 1, n * 2)
printf("%lld ", a[i][j]);
putchar(10);
}
4. 特征值与特征向量
有时，对一个向量进行一次线性变换（左乘一个矩阵），这个向量方向不变，只是长度发生了变化，即 (Av=lambda v)。也可以表示成 ((A-lambda I)v=0)
对于一个矩阵 (A)，若存在数 (lambda) 和向量 (v) 使得 (Av=lambda v)，称 (lambda) 为矩阵 (A) 的 特征值，(v) 为 特征向量。同时定义 (f(x)=|xI-A|) 为 特征多项式，(|xI-A|=0) 为 特征方程。
定理 (n) 阶矩阵 (A) 的 (n) 个特征值，为其特征方程的 (n) 个根。特征值 (lambda) 对应的特征向量为 ((lambda I-A)x=0) 的非零解。
注意一个特征值对应无穷多个特征向量，因为可以对向量进行数乘，仍满足定义。
矩阵对角化
相似矩阵：对于两个矩阵 (A,B)，若存在可逆矩阵 (P) 使得 (P^{-1}AP=B)，则 (A,B) 相似。
考虑矩阵 (A) 和对角矩阵 (Lambda=left[egin{matrix}lambda_1\&lambda_2\&&ddots\&&&lambda_nend{matrix} ight]) 相似的条件。设矩阵 (P) 由 (n) 个列向量构成即 (P=[p_1,p_2,dots,p_n])，且 (P) 可逆
则 (A[p_1,p_2,dots,p_n]=[p_1,p_2,dots,p_n]Lambda) 即 ([Ap_1,Ap_2,dots,Ap_n]=[lambda_1p_1,lambda_2p_2,dots,lambda_np_n])。
那么 (Ap_i=lambda_ip_i)。这不就是特征值和特征向量的定义吗？
于是有 (A=PLambda P^{-1})，其中 (Lambda) 是由 (A) 的所有特征值构成的对角矩阵，(P) 的第 (i) 列是特征值 (lambda_i) 对应的特征向量。
于是我们可以快速计算一些矩阵的幂，(A^n=(PLambda P^{-1})^n=PLambda^nP^{-1})。
递推数列通项公式
这里只讨论二阶齐次线性递推数列。
首先我们可以直接用上文的结论。比如斐波纳契转移矩阵 (left[egin{matrix}1~1\1~0end{matrix}
ight])，特征方程 (left|egin{matrix}lambda-1~-1\-1~lambdaend{matrix}
ight|=0) 即 (lambda^2-lambda-1=0)，解得 (lambda=dfrac{1pmsqrt 5}{2})，并分别求出特征向量 ([1,dfrac{sqrt 5-1}{2}]^T,[1,-dfrac{sqrt 5+1}{2}]^T)。然后一顿爆算求出 (P) 的逆 (left[egin{matrix}dfrac{5+sqrt5}{10}~dfrac{sqrt5}{5}\dfrac{5-sqrt5}{10}~-dfrac{sqrt5}{5}end{matrix}
ight])。把 (P,Lambda^n,P^{-1}) 依次乘起来，最后乘上 ([1,0]^T)，就可得出通项公式 (f_n=dfrac{sqrt 5}{5}left[(dfrac{1+sqrt5}{2})^n-(dfrac{1-sqrt5}{2})^n
ight])。
这么麻烦，为什么不用生成函数呢
考虑一个数列 (f_i=pf_{i-1}+qf_{i-2})，其中 (f_0,f_1) 已知
尝试构造等比数列：设存在 (a,b) 使得 (f_{i+1}-af_i=b(f_i-af_{i-1}))，即 (f_{i+1}=(a+b)f_i-abf_{i-1})
发现 (a,b) 是方程 (x^2-px-q=0) 的根。而这个方程就是转移矩阵的特征方程
所以 ({f_{i+1}-af_i}) 是公比为 (b) 的等比数列，有 (f_{i+1}-af_i=b^i(f_1-af_0))。
又发现 (a,b) 也满足 (f_{i+1}-bf_i=a(f_i-bf_{i-1}))
所以 ({f_{i+1}-bf_i}) 是公比为 (a) 的等差数列，有 (f_{i+1}-bf_i=a^i(f_1-bf_0))。
两式相减，得到 (f_i=dfrac{f_1-af_0}{b-a}b^i+dfrac{f_1-bf_0}{a-b}a^i)
其实并不用每次推。可以设 (f_i=alpha a^i+eta b^i)，解一个方程就完了
比如斐波纳契数列，(a=dfrac{1+sqrt 5}{2},b=dfrac{1-sqrt 5}{2})，设 (f_i=alpha a^i+eta b^i)，则 (egin{cases}f_0=alpha +eta=0\f_1=frac{1+sqrt 5}{2}alpha+frac{1-sqrt 5}{2}etaend{cases})
解得 (alpha=dfrac{sqrt 5}{5},eta=-dfrac{sqrt 5}{5})，那么 (f_n=dfrac{sqrt 5}{5}left[(dfrac{1+sqrt5}{2})^n-(dfrac{1-sqrt5}{2})^n ight])。
对于更高阶的情况，记个结论：设 (f_i=alpha a^i+eta b^i+gamma c^i+dots)，列方程组
常系数齐次线性递推
首先可以写出一个 (k imes k) 的转移矩阵：(A=left[egin{matrix}f_1~f_2~f_3~dots f_{k-1}~f_k\1~~~0~~~0~~dots~~~0~~~0\0~~~1~~~0~~dots~~~0~~~0\0~~~0~~~1~~dots~~~0~~~0\vdots~~~vdots~~~vdots~~~ddots~~~vdots~~~vdots\0~~~0~~~0~~dots~~~1~~~0end{matrix} ight])
初始矩阵 (B=[a_{k-1},a_{k-2},dots,a_0]^T)，答案是 ((A^nB)_k)。
凯莱-哈密顿定理：对于矩阵 (A)，(f(x)=|xI-A|) 为其特征多项式，有 (f(A)=O)，其中 (O) 是零矩阵。
本题中，可以用行列式性质 7 和 8 计算出 (f(x)=x^k-f_1x^{k-1}-f_2^{k-2}-dots-f_k)
那么只需计算 (g(x)=x^nmod f(x)=sumlimits_{i=0}^{k-1}g_ix^i)，答案为 ((g(A)B)_k=(sumlimits_{i=0}^{k-1}g_iA^iB)_k=sumlimits_{i=0}^{k-1}g_i(A^iB)_k=sumlimits_{i=0}^{k-1}g_ia_i)。
计算 (g(x)) 可以用倍增 + 多项式取模。时间复杂度 (O(klog klog n))。