数论 矩阵

矩阵

定义

由 \(m\times n\) 个数 \(a_{ij}\) 排成的 \(m\) 行 \(n\) 列的数表称为 \(m\) 行 \(n\) 列的矩阵，简称 \(m\times n\) 矩阵。记作：

\[A=\left[\begin{matrix} a_{11}&a_{12}&\dots&a_{1n}\\ a_{21}&a_{22}&\dots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}&a_{m2}&\dots&a_{mn} \end{matrix}\right] \]

这 \(m\times n\) 个数称为矩阵 \(A\) 的元素，简称元，数 \(a_{ij}\) 位于矩阵 \(A\) 的第 \(i\) 行第 \(j\) 列

这个矩阵还可以简记为： \(A=A_{m\times n}=\left(a_{ij}\right)_{m\times n}=\left(a_{ij}\right)\)

特殊矩阵

• 行列均为 \(n\) 的矩阵称为 \(n\) 阶方阵（也叫 \(n\) 阶矩阵）

• 只有一行的矩阵称为行矩阵（行向量）

• 只有一列的矩阵称为列矩阵（列向量）

• 若两个矩阵的行数和列数分别相同，则称这两个矩阵为同型矩阵

• 元素全为 \(0\) 的矩阵称为零矩阵，记作 \(0\) 。

  • 不同阶数的零矩阵不相同

• 除主对角线外的元素均为 \(0\) 的矩阵称为对角矩阵

  \[A=\left[ \begin{matrix} \lambda_1&0&\cdots&0\\ 0&\lambda_2&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&\lambda_n \end{matrix} \right]=diag(\lambda_1,\lambda_2,\dots,\lambda_n) \]

• 主对角线的元素为 \(1\) ，其他元素都为 \(0\) 的 \(n\) 阶矩阵称为单位矩阵

• 若两个同型矩阵的对应元素相等，则两个矩阵相等

运算

矩阵加减

• 只有两个同型矩阵才可以做加减法

\[设\ A=\left(a_{ij}\right),B=\left(b_{ij}\right)\\ A\pm B=\left[ \begin{matrix} a_{11}\pm b_{11}&a_{12}\pm b_{12}&\cdots&a_{1n}\pm b_{1n}\\ a_{21}\pm b_{21}&a_{22}\pm b_{22}&\cdots&a_{2n}\pm b_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}\pm b_{m1}&a_{m2}\pm b_{m2}&\cdots&a_{mn}\pm b_{mn} \end{matrix} \right] \]

• 性质
  • 交换律： \(A+B=B+A\)
  • 结合律： \((A+B)+C=A+(B+C)\)

矩阵乘法

• 设 $A_{m\times n} ,B_{a\times b} $ ，当且仅当 \(n=a\) 时，两个矩阵可以以 \(A\times B\) 的形式相乘

\[若 A_{m\times n}=(a_{ij}),B_{a\times b}={b_{ij}}\\ 设 C_{m\times b}=A\times B=(c_{ij})\\ 有 c_{ij}=a_{i1}b_{1j}+a_{i2}b_{2j}+\cdots+a_{in}b_{aj} \]

• 性质
  • 不满足交换律
  • 结合律： \((AB)C=A(BC)\)
  • 左分配律： \((A+B)C=AC+BC\)
  • 右分配律： \(C(A+B)=CA+CB\)
  • 数乘： \(k(AB)=(kA)B=A(kB)\)

代码实现

int a[maxn][maxn], b[maxn][maxn], c[maxn][maxn];
int main(){
    int m,n,k;
    scanf("%d%d%d", &m, &n, &k);
    for(int i = 1; i <= m; ++i)
        for(int j = 1; j <= n; ++j)
            scanf("%d", &a[i][j]);
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= k; ++j)
            scanf("%d", &b[i][j]);
    for(int i = 1; i <= m; ++i)//枚举a的行数m
        for(int j = 1; j <= k; ++j)//枚举b的列数k
            for(int kk = 1; kk <= n; ++kk)//枚举a的列(b的行)
                c[i][j] += a[i][kk]*b[kk][j];
    return 0;
}

高斯消元

定义

以下内容引自wikipedia

高斯消元法(Gaussian Elimination)是数学上线性代数中的一个算法，可以把矩阵转化为行阶梯形矩阵。高斯消元法可用来为线性方程组求解，求出矩阵的秩，以及求出可逆方阵的逆矩阵

该方法以数学家卡尔·高斯命名，但最早出现于《九章算术》

行阶梯形矩阵：满足下列条件的矩阵被称为行阶梯形矩阵：

• 元素不全为零的行（非零行）在所有元素全为零的行（全零行）的上面，即全零行都在矩阵的底部
• 非零行的首项系数，也称作主元，即最左边的首个非零元素，严格地比上面行的主元靠右
  • （前两条的推论）首项系数所在的列，在该首项系数下面的元素都是零

例：下面的这个 \(3\times 4\) 矩阵是行阶梯形矩阵：

\[\left[ {\begin{array}{c:c} \begin{matrix} 1&a_1&a_2\\ 0&2&a_3\\ 0&0&1 \end{matrix}& \begin{matrix} b_1\\b_2\\b_3 \end{matrix} \end{array}} \right] \]

矩阵的秩：在线性代数中，一个矩阵 \(A\) 的列秩是 \(A\) 的线性无关的纵列的极大数目。类似地，行秩是 \(A\) 的线性无关的横行的极大数目。同一行秩和列秩总是相等的，因此可以简单地称为矩阵 \(A\) 的秩，写作 \(\text{r}(A),\rank(A)或\text{rk}(A)\)

例：在下面这个 \(4\times4\) 矩阵中：

\[A=\left[ \begin{matrix} 2&4&1&3\\ -1&-2&1&0\\ 0&0&2&2\\ 3&6&2&5 \end{matrix} \right] \]

不难发现第二纵列是第一纵列的两倍，而第四纵列等于第一和第三纵列的总和。第一和第三纵列是线性无关的，所以 \(A\) 的秩是 \(2\) ，即 \(\rank(A)=2\)

• 性质：一个矩阵的行秩和列秩总是相等的（证明）

• 性质：矩阵的初等变换不改变矩阵的秩

逆矩阵：（又称乘法反方阵、反矩阵）在线性代数中，给定一个 \(n\) 阶方阵 \(A\) ，若存在一个 \(n\) 阶方阵 \(B\) ，使得 \(A\times B=B\times A =I_n\) ，其中 \(I_n\) 是 \(n\) 阶单位矩阵，则称 \(A\) 是可逆的，且 \(B\) 是 \(A\) 的逆矩阵，记作 \(A^{-1}=B\)

• 只有方阵（ \(n\times n\) 的矩阵）才可能有逆矩阵
• 若 \(A\) 的逆矩阵存在，则称 \(A\) 为非奇异方阵或可逆方阵

补充概念：

• 系数矩阵、增广矩阵：

  设有一个方程组 \(S\) ：

  \[\begin{equation} S=\left\{ \begin{aligned} a_{11}x_1+a_{12}x_2+\dots+a_{1n}&=b_1\\ a_{21}x_1+a_{22}x_2+\dots+a_{2n}&=b_2\\ &\vdots\\ a_{m1}x_1+a_{m2}x_2+\dots+a_{mn}&=b_m \end{aligned} \right. \end{equation} \]

  那么下面的这个 \(n\times m\) 的矩阵称为上面的方程组的系数矩阵：

  \[A=\left[ \begin{matrix} a_{11}&a_{12}&\dots&a_{1n}\\ a_{21}&a_{22}&\dots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}&a_{m2}&\dots&a_{mn} \end{matrix} \right] \]

  而增广矩阵就是在系数矩阵的右边添上一列，这一列是线性方程组的等号右边的值：

  \[B=\left[ {\begin{array}{c:c} \begin{matrix} a_{11}&a_{12}&\dots&a_{1n}\\ a_{21}&a_{22}&\dots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}&a_{m2}&\dots&a_{mn} \end{matrix}& \begin{matrix} b_1\\b_2\\\vdots\\b_m \end{matrix} \end{array}} \right] \]

基本原理

已知线性方程组：

\[\begin{equation} S=\left\{ \begin{aligned} 2x+y-z&=8\\ -3x-y+2z&=-11\\ -2x+y+2z&=-3 \end{aligned} \right. \end{equation} \]

这个方程组可以变换为如下形式：

\[\begin{equation} S=\left\{ \begin{aligned} x+\frac12y-\frac12z&=4\\ -x-\frac13y+\frac23z&=-\frac{11}3\\ -x+\frac12y+z&=-\frac32 \end{aligned} \right. \end{equation} \]

\[\Rightarrow\begin{equation} S=\left\{ \begin{aligned} x+\frac12y-\frac12z&=4\\ y+z&=2\\ -y-\frac12z&=-\frac52 \end{aligned} \right. \end{equation} \]

\[\Rightarrow\begin{equation} S=\left\{ \begin{aligned} x+\frac12y-\frac12z&=4\\ y+z&=2\\ z&=-1 \end{aligned} \right. \end{equation} \]

再将 \(z\) 回代，我们可以依次解出 \(y,x\) 的值

在上述过程中，方程组 \(S\) 的增广矩阵变化如下：

\[\left[ \begin{array}{c:c} \begin{matrix} 2&1&-1\\ -3&-1&2\\ -2&1&2 \end{matrix}& \begin{matrix} 8\\-11\\-3 \end{matrix} \end{array} \right]\\ \Rightarrow\left[ \begin{array}{c:c} \begin{matrix} 1&\frac12&-\frac12\\ -1&-\frac13&\frac23\\ -1&\frac12&1 \end{matrix}& \begin{matrix} 4\\-\frac{11}3\\-\frac32 \end{matrix} \end{array} \right]\\ \Rightarrow\left[ \begin{array}{c:c} \begin{matrix} 1&\frac12&-\frac12\\ 0&1&1\\ 0&-1&-\frac12 \end{matrix}& \begin{matrix} 4\\2\\-\frac52 \end{matrix} \end{array} \right]\\ \Rightarrow\left[ \begin{array}{c:c} \begin{matrix} 1&\frac12&-\frac12\\ -1&-\frac13&\frac23\\ -1&\frac12&1 \end{matrix}& \begin{matrix} 4\\-\frac{11}3\\-\frac32 \end{matrix} \end{array} \right]\\ \Rightarrow\left[ \begin{array}{c:c} \begin{matrix} 1&\frac12&-\frac12\\ 0&1&1\\ 0&0&1 \end{matrix}& \begin{matrix} 4\\2\\-1 \end{matrix} \end{array} \right] \]

我们可以发现，这个过程实际上就是矩阵的初等变换，将一个增广矩阵变换为了一个行阶梯形矩阵。之后再将已知的答案一个个地带入已经被简化的等式中的未知数中，就可以求出其余的答案了。

方程组解的判断

无解

当矩阵中出现某一行形如 \((0,0,\dots,a),a\not=0\) ，即此时有 \(0=a,a\not=0\) ，显然矛盾，方程组无解

唯一解

当行阶梯形矩阵形成严格的上（下）三角阵，或者可以化为 \(n\) 阶最简矩阵，那么有唯一解

注： \(n\) 阶最简矩阵（简化行阶梯形矩阵）：用上面的例子，上面的矩阵的 \(n\) 阶最简矩阵为：

\[\left[ \begin{array}{c:c} \begin{matrix} 1&0&0\\0&1&0\\0&0&1 \end{matrix}& \begin{matrix} 2\\3\\-1 \end{matrix} \end{array} \right] \]

无穷解

当矩阵中出现至少一行形如 \((0,0,\dots,0,a),a=0\) 时，有无穷解

每出现一个全零行就有一个自由元，有 \(k\) 个全零行就有 \(k\) 个自由元

代码

e.g. 高斯消元（模板）

描述

已知 \(n\) 元线性一次方程组 \(S\) ：

\[\begin{equation} S=\left\{ \begin{aligned} a_{11}x_1+a_{12}x_2+\dots+a_{1n}x_n&=b_1\\ a_{21}x_1+a_{22}x_2+\dots+a_{2n}x_n&=b_2\\ \vdots&\\ a_{n1}x_1+a_{n2}x_2+\dots+a_{nn}x_n&=b_n\\ \end{aligned} \right. \end{equation} \]

求 \(S\) 的解，若不存在唯一解输出 No Solution

I/O

Input: 第一行一个正整数 \(n\)

之后 \(n\) 行每行 \(n+1\) 个整数，分别代表 \(a_{i1},a_{i2},\dots,a_{in},b_i\)

Output: 若存在唯一解，输出共 \(n\) 行，每行一个数 \(x_i\) （保留两位小数）

若不存在唯一解输出一行 No Solution

范围

\(1\le n\le100,\ |a_{ij}|,|b_i|\lt 10^4\)

解法

\(O(n^3)\)

#include<bits/stdc++.h>
using namespace std;
#define GO(u,v,i) for(register int i=u;i<=v;i++)
template<class t>inline t fr(){
    register t num=0,dis=1;
    register char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')dis=-1;ch=getchar();}
    while(ch<='9'&&ch>='0'){num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
    return num*dis;
}
template<class t>inline void fw(t num){
    if(num>9)fw(num/10);
    putchar(num%10+'0');
}
template<class t>inline void fw(t num,char ch){
    if(num<0)num=-num,putchar('-');
    fw(num);putchar(ch);
}
typedef double lf;
const int maxn=100+12;
const lf eps=1e-6;//误差
lf a[maxn][maxn];//增广矩阵
int n;//题目条件
//推荐阅读顺序：main()->gauss()->prnt()
inline bool gauss(lf a[][maxn],int n){
    for(int i=1,maxi;i<=n;i++){//依次遍历n个方程
        maxi=i;//记录第i列i+1行到n行间系数绝对值最大行
        GO(i+1,n,j)if(fabs(a[j][i])>fabs(a[maxi][i]))maxi=j;
        //若第i~n行的第i列均为0，说明方程组有无穷解或无解
        if(fabs(a[maxi][i])<eps)return false;
        if(maxi^i)GO(1,n+1,j)swap(a[i][j],a[maxi][j]);//交换两行
        //消元
        GO(i+1,n,j){
            lf tmp=a[j][i]/a[i][i];
            if(fabs(a[j][i])>eps)GO(i,n+1,k)a[j][k]-=tmp*a[i][k];
        }
    return true;
}

inline void prnt(lf a[][maxn],int n){
    //求解输出
    for(int i=n;i>0;i--){//回代求解
        GO(i+1,n,j)a[i][n+1]-=a[j][n+1]*a[i][j];//a[i][i]已化为1
        a[i][n+1]/=a[i][i];
        if(fabs(a[i][n+1])<eps)a[i][n+1]=0;//防负
    }
    GO(1,n,i)printf("%.2lf\n",a[i][n+1]);
}

signed main(){
    n=fr<int>();
    GO(1,n,i)GO(1,n+1,j)scanf("%lf",&a[i][j]);
    if(gauss(a,n))prnt(a,n);
    else printf("No Solution\n");
    return 0;
}

高斯-约旦消元法

高斯-约旦消元法（亦称高斯-若尔当消元法）在高斯消元法的基础上，将有唯一解的矩阵进一步约简为简化行阶梯形矩阵，因而无需回代求解。

#include<bits/stdc++.h>
using namespace std;
#define GO(u,v,i) for(register int i=u;i<=v;i++)
template<class t>inline t fr(){
    register t num=0,dis=1;
    register char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')dis=-1;ch=getchar();}
    while(ch<='9'&&ch>='0'){num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
    return num*dis;
}
template<class t>inline void fw(t num){
    if(num>9)fw(num/10);
    putchar(num%10+'0');
}
template<class t>inline void fw(t num,char ch){
    if(num<0)num=-num,putchar('-');
    fw(num);putchar(ch);
}
typedef double lf;
const int maxn=100+12;
const lf eps=1e-6;//误差
lf a[maxn][maxn];//增广矩阵
int n;//题目条件
//推荐阅读顺序：main()->gauss()->prnt()
inline bool gauss(lf a[][maxn],int n){
    for(int i=1,maxi;i<=n;i++){
        //枚举每个方程
        maxi=i;
        GO(i+1,n,j)if(fabs(a[j][i])>fabs(a[maxi][i]))maxi=j;
        if(fabs(a[maxi][i])<eps)return false;//没有唯一解
        if(i^maxi)//fabs(a[i][i])非最大，交换
            GO(i,n+1,j)swap(a[i][j],a[maxi][j]);
        lf tmp=a[i][i];//抽出系数，否则变成1以后就没法变了
        GO(i,n+1,j)a[i][j]/=tmp;//第i行除以a[i][i]，变成1
        GO(1,n,j){
            if(i!=j){
            	tmp=a[j][i];
            	GO(i,n+1,k)a[j][k]-=tmp*a[i][k];
        	}
        }
    }
    return true;
}

inline void prnt(lf a[][maxn],int n){
    //输出
    GO(1,n,i)printf("%.2lf\n",a[i][n+1]);
}

signed main(){
    n=fr<int>();
    GO(1,n,i)GO(1,n+1,j)scanf("%lf",&a[i][j]);
    if(gauss(a,n))prnt(a,n);
    else printf("No Solution\n");
    return 0;
}

高斯-约旦消元法的常数稍大

e.g. 球形空间产生器

题目传送门：P4035 （省选/NOI-）

这道题的难点在于如何将这个连等式转化为方程

我们假设球心的坐标为 \((x_1,x_2,\dots,x_n)\) ，给出的点的坐标依次为 \((a_1,a_2,\dots,a_n),(b_1,b_2,\dots,b_n),\dots,(k_1,k_2,\dots,k_n)\)

根据题目提示，很容易想到以下等式：

\[(a_1-x_1)^2+(a_2-x_2)^2+\dots+(a_n-x_n)^2\\ =(b_1-x_1)^2+(b_2-x_2)^2+\dots+(b_n-x_n)^2\\ =\dots=(k_1-x_1)^2+\dots+(k_n-x_n)^2 \]

我们拆分出其中的两个多项式：

\[(a_1-x_1)^2+(a_2-x_2)^2+\dots+(a_n-x_n)^2\\ =(b_1-x_1)^2+(b_2-x_2)^2+\dots+(b_n-x_n)^2 \]

开括号

\[a_1^2-2a_1x_1+x_1^2+a_2^2-2a_2x_2+x_2^2+\dots+a_n^2-2a_nx_n+x_n^2\\ =b_1^2-2b_1x_1+x_1^2+b_2^2-2b_2x_2+x_2^2+\dots+b_n^2-2b_nx_n+x_n^2 \]

观察到等式两边都有\((x_1^2+x_2^2+\dots+x_n^2)\) ，消去，同时将常数项移至等式右边，将未知数移至等式左边

\[(2b_1-2a_1)x_1+(2b_2-2a_2)x_2+\dots+(2b_n-2a_n)x_n\\ =a_1^2+b_1^2+a_2^2+b_2^2+\dots+a_n^2+b_n^2 \]

可以发现，对于这 \((n+1)\) 个连等的多项式，我们可以每次取出其中的相邻的两个多项式，然后仿照上面的方法化简，最终得到 \(n\) 个等式构成的 \(n\) 元一次方程组

而解这个方程组则可以使用高斯消元法

#include<bits/stdc++.h>
using namespace std;
#define GO(u,v,i) for(register int i=u;i<=v;i++)
template<class t>inline t fr(){
    register t num=0,dis=1;
    register char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')dis=-1;ch=getchar();}
    while(ch<='9'&&ch>='0'){num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
    return num*dis;
}
template<class t>inline void fw(t num){
    if(num>9)fw(num/10);
    putchar(num%10+'0');
}
template<class t>inline void fw(t num,char ch){
    if(num<0)num=-num,putchar('-');
    fw(num);putchar(ch);
}
typedef double lf;
const int maxn=10+12;
const lf eps=1e-6;
int n;
lf a[maxn][maxn];

#include<bits/stdc++.h>
using namespace std;
#define GO(u,v,i) for(register int i=u;i<=v;i++)
template<class t>inline t fr(){
    register t num=0,dis=1;
    register char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')dis=-1;ch=getchar();}
    while(ch<='9'&&ch>='0'){num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
    return num*dis;
}
template<class t>inline void fw(t num){
    if(num>9)fw(num/10);
    putchar(num%10+'0');
}
template<class t>inline void fw(t num,char ch){
    if(num<0)num=-num,putchar('-');
    fw(num);putchar(ch);
}
typedef double lf;
const int maxn=100+12;
const lf eps=1e-6;//误差
lf a[maxn][maxn];//增广矩阵
int n;//题目条件
//推荐阅读顺序：main()->gauss()->prnt()
inline bool gauss(lf a[][maxn],int n){
    for(int i=1,maxi;i<=n;i++){
        //枚举每个方程
        maxi=i;
        GO(i+1,n,j)if(fabs(a[j][i])>fabs(a[maxi][i]))maxi=j;
        if(fabs(a[maxi][i])<eps)return false;//没有唯一解
        if(i^maxi)//fabs(a[i][i])非最大，交换
            GO(i,n+1,j)swap(a[i][j],a[maxi][j]);
        lf tmp=a[i][i];//抽出系数，否则变成1以后就没法变了
        GO(i,n+1,j)a[i][j]/=tmp;//第i行除以a[i][i]，变成1
        GO(1,n,j){
            if(i!=j){
            	tmp=a[j][i];
            	GO(i,n+1,k)a[j][k]-=tmp*a[i][k];
        	}
        }
    }
    return true;
}

inline void prnt(lf a[][maxn],int n){
    //输出
    GO(1,n,i)printf("%.2lf\n",a[i][n+1]);
}

signed main(){
    n=fr<int>();
    GO(1,n,i)GO(1,n+1,j)scanf("%lf",&a[i][j]);
    if(gauss(a,n))prnt(a,n);
    else printf("No Solution\n");
    return 0;
}

signed main(){
    n=fr<int>();
    GO(1,n+1,i)GO(1,n,j)scanf("%lf",&a[i][j]);
    GO(1,n,i){
        GO(1,n,j){//常数项
            a[i][n+1]-=a[i][j]*a[i][j];
        	a[i][n+1]+=a[i+1][j]*a[i+1][j];
        }
        GO(1,n,j){//各项系数
            a[i][j]=-2*a[i][j];
        	a[i][j]+=2*a[i+1][j];
        }
    }
    gauss(a,n);//高斯消元
    prnt(a,n);//求解输出
    return 0;
}

判断无解与无穷解

常规的高斯消元与高斯-约旦消元在遇到无解和无穷解时都会直接结束，无法进行判断。那要如何判断无解和无穷解呢？

根据上文推导，在把增广矩阵消解成上三角型的行阶梯形矩阵后，若存在 \(a_{i,i}=0,a_{i,n+1}=0\) 时有无穷解，若存在一行 \(a_{i,i}=0,a_{i,n+1}\not=0\) 时无解，但用正常的消解法无法做出准确判断，例如下面的矩阵：

\[S=\left[ \begin{array}{c:c} \begin{matrix} 1&2&3&4&5&6\\ 0&0&1&2&3&4\\ 0&0&0&2&6&7\\ 0&0&0&0&6&8\\ 0&0&0&0&0&1\\ 0&0&0&0&0&0 \end{matrix}& \begin{matrix} 10\\15\\30\\10\\30\\0 \end{matrix} \end{array} \right] \]

这是一个行阶梯形矩阵，其中 \(a_{5,5}=0\) ，但 \(a_{5,7}=30\) ，判断方程无解，但实际上是无穷解

我们交换这个增广矩阵的第一列和第六列，这样不会影响方程的解

\[S=\left[ \begin{array}{c:c} \begin{matrix} 6&2&3&4&5&1\\ 4&0&1&2&3&0\\ 7&0&0&2&6&0\\ 8&0&0&0&6&0\\ 1&0&0&0&0&0\\ 0&0&0&0&0&0 \end{matrix}& \begin{matrix} 10\\15\\30\\10\\30\\0 \end{matrix} \end{array} \right] \]

这个矩阵可以消解为下面这个矩阵（为便于观察，对实数取整）：

\[S=\left[ \begin{array}{c:c} \begin{matrix} 8&0&0&0&6&0\\ 0&2&3&4&0&1\\ 0&0&1&2&0&0\\ 0&0&0&2&1&0\\ 0&0&0&0&-1&0\\ 0&0&0&0&0&0 \end{matrix}& \begin{matrix} 10\\2\\10\\21\\29\\0 \end{matrix} \end{array} \right] \]

这个矩阵除了 \(a_{6,6}=0\) 外其他 \(a_{i,i}\not=0\) ，所以方程组无穷解

根据上面的例子可以看出，方程的顺序会影响解的判断。那如何在不改变变量位置的情况下能正确地对防方程组进行判断呢？

我们对高斯消元法稍作改进

在消解时，把坡度一致的倒三角改成坡度不一致的倒三角。具体做法是：在消解第 \(i\) 行时，若 \(i\to n\) 行的第 \(i\) 列全为零，常规做法是枚举第 \(i+1\) 行消解 \(i+1\) 列，此时我们改变策略，我们保持消解第 \(i\) 行，但是消解第 \(i+1\) 列，若此时任找不到非 \(0\) 主元，就在第 \(i\) 行 \(i+2\) 列寻找非零主元，直到找到非零主元为止

e.g. 线性方程组

题目描述

已知 \(n\) 元一次方程组

\[\begin{equation} S=\left\{ \begin{aligned} a_{11}x_1+a_{12}x_2+\dots+a_{1n}x_n&=b_1\\ a_{21}x_1+a_{22}x_2+\dots+a_{2n}x_n&=b_2\\ \vdots&\\ a_{n1}x_1+a_{n2}x_2+\dots+a_{nn}x_n&=b_n\\ \end{aligned} \right. \end{equation} \]

若有唯一解，求解；若无解，输出 \(-1\) ；若有无穷解，输出 \(0\)

I/O

Input: 第一行一个正整数 \(n\) 表示方程数

接下来 \(n\) 行，每行 \(n+1\) 个整数，表示 \(a_{i1},a_{i2},\dots,a_{in},b_i\)

Output: 若无解，输出 \(-1\) ；若有无穷解，输出 \(0\)

若有唯一解，按如下格式输出（保留两位有效数字）

x1=1.00
x2=0.00
...
xn=8.21

范围

\(1\le n\le 50,|a_{ij}|\le 100,0\lt b_i\lt 300\)

解法

\(O(n^3)\)

#include<bits/stdc++.h>
using namespace std;
#define GO(u,v,i) for(register int i=u;i<=v;i++)
template<class t>inline t fr(){
    register t num=0,dis=1;
    register char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')dis=-1;ch=getchar();}
    while(ch<='9'&&ch>='0'){num=(num<<1)+(num<<3)+(ch^48);ch=getchar();}
    return num*dis;
}
template<class t>inline void fw(t num){
    if(num>9)fw(num/10);
    putchar(num%10+'0');
}
template<class t>inline void fw(t num,char ch){
    if(num<0)num=-num,putchar('-');
    fw(num);putchar(ch);
}
typedef double lf;
const int maxn=100+12;
const lf eps=1e-6;
lf a[maxn][maxn];
int n;
//推荐阅读顺序：main()->gauss()->prnt()
inline void prnt(lf a[][maxn],int base,int n){
    if(base<=n){
        GO(base,n,i){
            if(fabs(a[i][n+1])>eps){
                fw(-1,'\n');
                return;
            }
        }
        fw(0,'\n');return;
    }
    GO(1,n,i){
        if(fabs(a[i][n+1])<eps)a[i][n+1]=0;
        printf("x%d=%.2lf\n",i,a[i][n+1]);
    }
}

inline void gauss(lf a[][maxn],int n){
    int base=1,maxi;//枚举作为基元的行
    GO(1,n,i){
        maxi=base;
        GO(base+1,n,j)if(fabs(a[j][i])>fabs(a[maxi][i]))maxi=j;
        //若a[maxi][i]=0表示base~n行的第i列全为零，此时继续在base行消解
        if(fabs(a[maxi][i])<eps)continue;
        if(base^maxi)GO(i,n+1,j)swap(a[base][j],a[maxi][j]);
        lf tmp=a[base][i];
        GO(i,n+1,j){
            a[base][j]/=tmp;
        }
        GO(1,n,j){
            if(j==base)continue;
            tmp=a[j][i];
            GO(i,n+1,k)a[j][k]-=tmp*a[base][k];
        }
        base++;
    }
    prnt(a,base,n);
}

signed main(){
    n=fr<int>();
    GO(1,n,i)GO(1,n+1,j)scanf("%lf",&a[i][j]);
    gauss(a,n);
    return 0;
}