线性代数
本文仅发布于此博客和作者的洛谷博客，不允许任何人以任何形式转载，无论是否标明出处及作者。
声明：由于这是一个OI里的线性代数讲解而不是数学的讲解，所以这篇文章只会出现和OI高度相关的内容。
矩阵基础
矩阵是一个按照长方阵列排列的数集。记矩阵\( A_{n,m}= \left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m}\\ a_{2,1} & a_{2,2} & \cdots & a_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{matrix}\right).\)
单位矩阵
定义单位矩阵\( I_{n,n}=\left(\begin{matrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 \end{matrix}\right).\)
即除主对角线上元素为1外，其他位置上元素均为0的方阵。
零矩阵
定义零矩阵\( O_{n,m}=\left(\begin{matrix} 0 & 0 & \cdots & 0\\ 0 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 0 \end{matrix}\right).\)
全零的矩阵。
转置矩阵
定义\(A_{n,m}= \left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,m}\\ a_{2,1} & a_{2,2} & \cdots & a_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,m} \end{matrix}\right)\)的转置矩阵为\(A^T= \left(\begin{matrix} a_{1,1} & a_{2,1} & \cdots & a_{m,1}\\ a_{1,2} & a_{2,2} & \cdots & a_{m,2}\\ \vdots & \vdots & \ddots & \vdots\\ a_{1,n} & a_{2,n} & \cdots & a_{m,n} \end{matrix}\right)\)
即对\(A\)关于主对角线翻转。
虽然零矩阵和转置矩阵没什么用
矩阵加法
两个矩阵相加时，在每一位上相加即可。\( A_{n,m}+B_{n,m}= \left(\begin{matrix} a_{1,1}+b_{1,1} & a_{1,2}+b_{1,2} & \cdots & a_{1,m}+b_{1,m}\\ a_{2,1}+b_{2,1} & a_{2,2}+b_{2,2} & \cdots & a_{2,m}+b_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n,1}+b_{n,1} & a_{n,2}+b_{n,2} & \cdots & a_{n,m}+b_{n,m} \end{matrix}\right).\)
使用代码可能更好理解。
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
c[i][j]=a[i][j]+b[i][j];
}
}
矩阵减法
同理。
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
c[i][j]=a[i][j]-b[i][j];
}
}
注意，行数和列数相同的矩阵才能做加减法。
矩阵数乘
矩阵的数乘和乘法不同。矩阵\(A\)和数\(k\)相乘，相当于\(A\)的每一位和\(k\)相乘。
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
c[i][j]=a[i][j]*k;
}
}
易得，上面的运算都满足交换律，结合律，复杂度均为\(O(n^2)\)
矩阵乘法
给出两个矩阵\(A_{n,p},B_{p,m}\)，相乘后得到一个矩阵\(C_{n,m},\)，满足\(\large C_{i,j}=\sum\limits_{k=1}^{p} A_{i,k}\times B_{k,j}\)，即“左行右列”。
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
c[i][j]=0;
for(int k=1;p<=p;k++){
c[i][j]+=a[i][k]*b[k][j];
}
}
}
复杂度\(O(n^3)\)
注意，矩阵乘法必须保证\(A\)的列数等于\(B\)的行数，而且矩阵乘法并不满足交换律。
矩阵快速幂
方阵\(A_{n,n}\)的\(k\)次幂记为\(A^k\)，\(A^k=\underbrace{A\times A\times\cdots\times A}_{k\text{个}}.\)
特殊地，\(A^0\)为单位矩阵\(I_{n,n}\)。
暴力做时，复杂度\(O(n^3k)\)，显然过高，尝试对其进行优化。
可以使用快速幂的方法。下面是标准的快速幂代码：
//a^k mod P
while(k){
if(k&1){
ans*=a;
ans%=P;
}
a*=a;
a%=P;
k>>=1;
}
本质上是一个把\(k\)次幂拆分成一些\(2\)的倍数次幂相乘。如果连这个都看不懂了的话建议还是花点时间进行一次复健。
我们对矩阵求幂进行完全一模一样的操作即可，总复杂度乘上一个矩阵乘法，\(O(n^3\log k)\)
矩阵快速幂和上面的东西可以一起集成在一个程序里：
#include<bits/stdc++.h>
#define int long long
using namespace std;
int P=1e9+7;
void debug(){//不要在意这些细节
cout<<"???what the hell"<<endl;
exit(0);
}
int mod(int k){//严谨的取模
return ((k%P)+P)%P;
}
struct matrix{
int n;
int m;
int v[105][205];
matrix(int t1,int t2){//初始化全0
n=t1;
m=t2;
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
v[i][j]=0;
}
}
}
};
matrix I(int n){//构造单位矩阵
matrix c(n,n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
c.v[i][j]=0;
if(i==j){
c.v[i][j]=1;
}
}
}
return c;
}
matrix O(int n,int m){//构造零矩阵
matrix c(n,m);
return c;
}
matrix T(matrix c){//构造转置矩阵
for(int i=1;i<=max(c.n,c.m);i++){
for(int j=1;j<=i;j++){
swap(c.v[i][j],c.v[j][i]);
}
}
swap(c.n,c.m);
return c;
}
matrix operator+(matrix a,matrix b){//矩阵加
if(a.n!=b.n||a.m!=b.m){
debug();
}
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]+b.v[i][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator-(matrix a,matrix b){//矩阵减
if(a.n!=b.n||a.m!=b.m){
debug();
}
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]-b.v[i][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator*(matrix a,int b){//矩阵数乘
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]*b;
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator*(matrix a,matrix b){//矩阵乘
if(a.m!=b.n){
debug();
}
matrix c(a.n,b.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=0;
for(int k=1;k<=a.m;k++){
c.v[i][j]+=a.v[i][k]*b.v[k][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
}
return c;
}
matrix pow(matrix a,int k){//矩阵快速幂
if(a.n!=a.m){
debug();
}
matrix c(a.n,a.m);
c=I(a.n);
while(k!=0){
if(k&1){
c=c*a;
}
a=a*a;
k>>=1;
}
return c;
}
void print(matrix c){
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
cout<<c.v[i][j]<<' ';
}
cout<<endl;
}
return;
}
signed main(){
return 0;
}
一个比较不是很全的矩阵全家桶就做好了。
这份代码先放放，处理一下矩阵加速。
矩阵加速
定义斐波那契数列\(F_i=F_{i-1}+F_{i-2}\)，其中\(F_1=F_2=1.\)尝试在快于\(O(n)\)的时间内计算出\(F_n \mod 10^9+7\)。
构造一个矩阵\(\left(\begin{matrix} F_i & F_{i+1} \end{matrix}\right)\)，可以发现矩阵变成\(\left(\begin{matrix} F_{i+1} & F_{i+2} \end{matrix}\right)\)即\(\left(\begin{matrix} F_{i+1} & F_i+F_{i+1} \end{matrix}\right)\)仅需乘上一个\(\left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)\)，对于每次变换皆如此。
因此\(\left(\begin{matrix} F_{k-1} & F_{k} \end{matrix}\right) = \left(\begin{matrix} F_1 & F_2 \end{matrix}\right)\times \left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)^{k-2}\)，而后可得\(F_k.\)
\(\left(\begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix}\right)^{k-2}\)使用矩阵快速幂，即可在\(O(\log k)\)的复杂度下完成，所以整个过程为\(O(\log k)\)。
模板题过程相似，快速幂的矩阵变为了\(3\times 3\)的而已。
另外，可以尝试想一下这道题。
给出\(n\)个操作，对于每个操作，输入两个数\(l,r\)，计算斐波那契数列上区间\([l,r]\)的数字和\(\mod\space 10^{16}+2137\)。斐波那契数列下标从\(1\)开始，依次为\(1,1,2,3\cdots\).
对于\(100\%\)的数据：保证\(1\leq n\leq 10^6,1\leq l\leq r\leq 10^9\) 。
题解：
设斐波那契数列为\(F\),对于每个查询\([l,r]\)，通过矩阵加速递推\(O(\log l)\)地找到\(F_l\)和\(F_{l+1}\)。设\(F_l=a,F_{l+1}=b,len=r-l+1.\)另设数列\(G,G_i=G_{i-1}+G_{i-2}+1.\)特殊地，\(G_1=0,G_2=1\)。仍然使用矩阵加速,在\(O(\log len)\)的时间内求\(F_{len},G_{len}\)，区间\([l,r]\)之和即为\(F_{len}\times a+G_{len}\times b.\)
对于每一次操作，复杂度为\(O(\log r)\)，总共为\(O(n\log r)\).
高斯消元
初等行变换
对于任意矩阵\(A\),可以进行三种初等行变换：
-
倍加：把某一行所有元素乘上某个数并加到一行。（原来那行不变）
-
倍乘：某一行乘上非零数\(k\)。
-
对换：两行元素互换。
注意到，对换可以由倍加操作完成。（不会的去看CSP初赛学怎么不用辅助空间交换元素）（虽然这个性质并没有什么用处，，，）
解线性方程组
我们把给出的方程组变成矩阵的形式：\(\left(\begin{matrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} & | & b_1\\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} & | & b_2\\ \vdots & \vdots & \ddots & \vdots & | & \vdots\\ a_{n,1} & a_{n,2} & \cdots & a_{n,n} & | & b_n\\ \end{matrix}\right).\)
通过初等行变换把左边的\(a\)方阵变换成单位矩阵，则\(b_i\)就是未知数的值。
具体实现方面，对于第\(i\)列，找到绝对值最小且非零的\(a_{i,j}\)（为了少掉精度），把这一行和第\(i\)行交换。然后用这一行消掉所有其他行的第\(i\)个元素。
如果找不到非零\(a_{i,j}\)，就无解。
代码：
#include<bits/stdc++.h>
#define int long long
using namespace std;
struct matrix{
int n;
int m;
double v[105][205];//别忘了改double
matrix(int t1,int t2){//初始化全0
n=t1;
m=t2;
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
v[i][j]=0;
}
}
}
};
signed main(){
int n;
cin>>n;
matrix c(n,n+1);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
cin>>c.v[i][j];
}
}
for(int j=1;j<=c.m-1;j++){//对每一列进行消元
double max_=-0x3f3f3f3f;
int maxpos=0;
for(int i=j;i<=c.n;i++){//找到j列最大的a[i][j]（注意：j行之前的行已经是单位矩阵的形式了，不能碰）
if(c.v[i][j]>max_){
max_=c.v[i][j];
maxpos=i;
}
}
for(int k=1;k<=c.m;k++){//交换maxpos行和j行
swap(c.v[j][k],c.v[maxpos][k]);
}
if(c.v[j][j]==0){//无解判定
cout<<"No Solution"<<endl;
return 0;
}
for(int i=1;i<=c.n;i++){//对其他行进行消元
if(i==j){//当前行跳过
continue;
}
double tmp=c.v[i][j]/c.v[j][j];//第j行全体乘上tmp再减到当前行上
for(int k=1;k<=c.m;k++){
c.v[i][k]-=c.v[j][k]*tmp;
}
}
double tmp=c.v[j][j];
for(int k=1;k<=c.m;k++){//a[j][j]化为1
c.v[j][k]/=tmp;
}
}
for(int i=1;i<=c.n;i++){//快乐输出
cout<<fixed<<setprecision(2)<<c.v[i][c.m]<<endl;
}
return 0;
}
矩阵求逆
定义方阵\(A\)的逆矩阵\(A^{-1}\)，使得\(AA^{-1}=I\)。
\(A\)的逆矩阵不一定存在。
以\(A=\left(\begin{matrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{matrix}\right)\)为例，逆矩阵的求法如下：
1.在\(A\)的右边补上一个单位矩阵。
\(C=\left(\begin{matrix} 2 & 1 & 1 & | & 1 & 0 & 0\\ 1 & 2 & 1 & | & 0 & 1 & 0\\ 1 & 1 & 2 & | & 0 & 0 & 1 \end{matrix}\right)\)
- 对左边进行高斯消元，消成单位矩阵。如果无法消元即代表逆矩阵不存在。
\(C'=\left(\begin{matrix} 1 & 0 & 0 & | & \frac34 & -\frac14 & -\frac14\\ 0 & 1 & 0 & | & -\frac14 & \frac34 & -\frac14\\ 0 & 0 & 1 & | & -\frac14 & -\frac14 & \frac34 \end{matrix}\right)\)
- 此时右半边就是逆矩阵了。
\(A^{-1}=\left(\begin{matrix} \frac34 & -\frac14 & -\frac14\\ -\frac14 & \frac34 & -\frac14\\ -\frac14 & -\frac14 & \frac34 \end{matrix}\right).\)
注意模板题是输出模意义下的逆矩阵，把除法全部换成逆元（费马小定理）即可。
行列式求值
方阵\(A\)的行列式记为\(|A|\)或\(\det(A)\)，它具有一个数值，值的计算遵循一个比较复杂的公式。行列式也可以进行高斯消元，并且值也会相应变化。
-
倍加值不变。
-
倍乘值乘k。
-
对换（自己和自己换不算）值反号。
另外，“上三角矩阵”（即主对角线一下均为0）的行列式的值为对角线乘积。
根据这个，我们可以把矩阵消成“上三角矩阵”（消成单位矩阵一步到位也可以的）计算值就可以了。
实现时，当无法消元时，答案为0.
奉上一份常数爆炸根本过不去的全家桶（
（其实还会有WA，因为这道题模数不是质数。正解是在消元的时候进行“辗转相减”而不是直接逆元）
#include<bits/stdc++.h>
#define int long long
using namespace std;
int P=1e9+7;
void debug(){
cout<<"???what the hell"<<endl;
exit(0);
}
int mod(int k){
return ((k%P)+P)%P;
}
int rev(int a){
int ans=1;
int k=P-2;
while(k!=0){
if(k&1){
ans*=a;
ans=mod(ans);
}
a*=a;
a=mod(a);
k>>=1;
}
return ans;
}
struct matrix{
int n;
int m;
int v[105][105];
matrix(int t1,int t2){
n=t1;
m=t2;
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
v[i][j]=0;
}
}
}
};
matrix I(int n){
matrix c(n,n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
c.v[i][j]=0;
if(i==j){
c.v[i][j]=1;
}
}
}
return c;
}
matrix O(int n,int m){
matrix c(n,m);
return c;
}
matrix T(matrix c){
for(int i=1;i<=max(c.n,c.m);i++){
for(int j=1;j<=i;j++){
swap(c.v[i][j],c.v[j][i]);
}
}
swap(c.n,c.m);
return c;
}
bool operator==(matrix a,matrix b){
if(a.n!=b.n||a.m!=b.m){
return false;
}
for(int i=1;i<=a.n;i++){
for(int j=1;j<=a.m;j++){
if(a.v[i][j]!=b.v[i][j]){
return false;
}
}
}
return true;
}
bool operator!=(matrix a,matrix b){
return !(a==b);
}
matrix operator+(matrix a,matrix b){
if(a.n!=b.n||a.m!=b.m){
debug();
}
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]+b.v[i][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator-(matrix a,matrix b){
if(a.n!=b.n||a.m!=b.m){
debug();
}
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]-b.v[i][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator*(matrix a,int b){
matrix c(a.n,a.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=a.v[i][j]*b;
c.v[i][j]=mod(c.v[i][j]);
}
}
return c;
}
matrix operator*(matrix a,matrix b){
if(a.m!=b.n){
debug();
}
matrix c(a.n,b.m);
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
c.v[i][j]=0;
for(int k=1;k<=a.m;k++){
c.v[i][j]+=a.v[i][k]*b.v[k][j];
c.v[i][j]=mod(c.v[i][j]);
}
}
}
return c;
}
matrix pow(matrix a,int k){
if(a.n!=a.m){
debug();
}
matrix c(a.n,a.m);
c=I(a.n);
while(k!=0){
if(k&1){
c=c*a;
}
a=a*a;
k>>=1;
}
return c;
}
pair<matrix,int> Gauss_jordan(matrix c){
int ans=1;//行列式求值专用
if(c.m<c.n){
debug();
}
for(int j=1;j<=c.n;j++){
int max_=-0x7fffffff;
int maxpos=0;
for(int i=j;i<=c.n;i++){
if(c.v[i][j]>max_&&c.v[i][j]!=0){
max_=c.v[i][j];
maxpos=i;
}
}
if(maxpos!=j){
ans*=-1;//换行就换号
ans=mod(ans);
}
for(int k=1;k<=c.m;k++){
swap(c.v[j][k],c.v[maxpos][k]);
}
if(c.v[j][j]==0){
matrix ERR(c.n,c.m);
return pair<matrix,int>{ERR,-0x7fffffff};
}
for(int i=1;i<=c.n;i++){
if(i==j){
continue;
}
int tmp=c.v[i][j]*rev(c.v[j][j]);
tmp=mod(tmp);
for(int k=1;k<=c.m;k++){
c.v[i][k]-=c.v[j][k]*tmp;
c.v[i][k]=mod(c.v[i][k]);
}
}
int tmp=rev(c.v[j][j]);
ans*=c.v[j][j];//行列式/k,答案*k
ans=mod(ans);
for(int k=1;k<=c.m;k++){
c.v[j][k]*=tmp;
c.v[j][k]=mod(c.v[j][k]);
}
}
return pair<matrix,int>{c,ans};
}
matrix rev(matrix a){
if(a.n!=a.m){
debug();
}
matrix c(a.n,a.m*2);//补单位矩阵
for(int i=1;i<=a.n;i++){
for(int j=1;j<=a.m;j++){
c.v[i][j]=a.v[i][j];
}
}
matrix tmp=I(a.n);
for(int i=1;i<=tmp.n;i++){
for(int j=1;j<=tmp.m;j++){
c.v[i][j+tmp.n]=tmp.v[i][j];
}
}
c=Gauss_jordan(c).first;//得到消元后矩阵
matrix ans(a.n,a.m);
for(int i=1;i<=a.n;i++){
for(int j=1;j<=a.m;j++){
a.v[i][j]=c.v[i][j+a.n];
}
}
return a;//拿到右半边
}
int det(matrix c){//直接集成在高斯消元里面算
return Gauss_jordan(c).second;
}
void scan(matrix &c){
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
cin>>c.v[i][j];
}
}
return;
}
void print(matrix c){
for(int i=1;i<=c.n;i++){
for(int j=1;j<=c.m;j++){
cout<<c.v[i][j]<<' ';
}
cout<<endl;
}
return;
}
signed main(){
return 0;
}
浙公网安备 33010602011771号