高斯消元法

高斯消元问题通常用来对线性方程组的求解

根据线性代数中所学，线性方程通常可表示成

a[i][0]*x0 + a[i][1]*x1 +a [i][2]*x2 ... + a[i][n-1]*xn-1 = b[i]

0<=i<n

这样就可以得到一个由n个线性方程组成方程组

我们通常将a[][]放入数组中的对应位置，b[] 作为一整列插入到a数组的最后一排

 

 

不断通过对行进行操作

比如 ri - k*rj , 用第i行的元素每一个都减去第j行对应列的元素的k倍，我们总是希望能够通过这个操作不断将矩阵朝着除去 b[]后的上三角

矩阵化简，就是令下三角矩阵均为0

那么就是不断从第一列开始不断将需要化简为0的项转化为0即可

那么我们就可以利用回代的方法，先算出最后一项x，然后不断往前递推，在第i行得到x[i]

 

 

当然这是最好的情况，解决线性方程组往往给的方程的数量和对应x的变量数量不一定相同

我令方程数量为equ ， 变量数量为 var ， 那么如果equ < var 的时候必然是有解的

因为至少有 var-equ+1个自由变量

因为每一行都会解决至多一个自由变量，当然如果某一行和前面的一行 的系数都存在倍数关系 exp: ai = k*aj

那么那一行其实效果是没有的，这个时候就一个自由变量都没有解决了

 

另外如果在化简矩阵的过程中某一行出现了（0,0,0,0....,0,b[i])（b[i] != 0)的情况那么就说明是不存在解的

 

所以对于我们来说，只要除去那些有相同效果的方程，剩下的方程组含有的方程为equ，那么

自由变量的个数就是 var - equ + 1 ， 当自由变量个数大于1的时候我们就可以认为是有无穷多个解的

 

我们利用这个思想解决计算机上的问题时，可以套用如下模版：

 

int gcd(int a , int b)
{
    if(b == 0) return a;
    else return gcd(b , a%b);
}

int lcm(int a , int b)
{
    return a/gcd(a , b)*b;//先除后乘防止溢出
}

inline int my_abs(int x)
{
    return x>=0?x:-x;
}

inline void my_swap(int &a , int &b)
{
    int tmp = a;
    a = b , b = tmp;
}
/*
高斯消元解线性方程组的解
返回-2表示有浮点数的解，但无整数解
返回-1表示无解
返回0表示唯一解
大于0，表示无穷解，返回的是自由变元个数
有equ个方程，var个变元。增广矩阵行数为equ，分别为0到equ-1，列数为var+1，分别为0到var
总是考虑上三角部分，所以内部更新总是只更新了上三角部分
*/
int gauss(int equ , int var)
{
    int max_r; // 当前这列绝对值最大的行
    int col; // 当前处理的列
    int ta , tb , i , j , k;
    int Lcm , tmp , free_x_num , free_index;

    for(i=0 ; i<=var ; i++)
        x[i] = 0 , free_x[i] = true;

    //转化为阶梯形矩阵
    col = 0;//一开始处理到第0列
    for(k=0 ; k<equ && col < var; k++,col++){
    /*
    枚举当前处理的行
    找到该行col列元素绝对值最大的那行与第k行交换，
    这样进行除法运算的时候就可以避免小数除以大数了
    */
        max_r = k;
        for(i = k+1 ; i<equ ; i++)
            if(my_abs(a[i][col]) > my_abs(a[max_r][col]))
                max_r = i;
        if(max_r != k){
            //交换两行的数据
            for(i=k ; i<var+1 ; i++)
                my_swap(a[max_r][i] , a[k][i]);
        }

        if(a[k][col] == 0){
            //说明该col列第k行以下全是0了，则处理当前行的下一列
            k--;
            continue;
        }

        for(i=k+1 ; i<equ ; i++){
            //枚举要删去的行
            if(a[i][col] != 0 ){
                Lcm = lcm(my_abs(a[i][col]) , my_abs(a[k][col]));
                ta = Lcm / my_abs(a[i][col]);
                tb = Lcm / my_abs(a[k][col]);
                if(a[i][col] * a[k][col] < 0) tb = -tb;//异号相加
                for(j = col ; j<var+1 ; j++)
                    a[i][j] = a[i][j]*ta - a[k][j]*tb ;
            }
        }
    }

    //1.无解的情况
    for(i=k ; i<equ ; i++){
        if(a[i][col] != 0) return -1;
    }

    //2.无穷解：在var*(var+1) 的增广阵中出现了(0,0,0....,0)这样的行，也即说明没有出现严格的上三角
    //出现这样情况的行数也就是自由变元的个数
    if(k < var){
        //首先初始化自由变元有var-k个，即不确定的变元至少有var-k个
        for(i=k-1 ; i>=0 ; i--){
            //第i行一定不会是(0,0,0,...,0)的情况，因为这样的行实在第k行到equ行
            //同样，第 i 行一定不会是(0,0 ... , a) , a != 0的情况

            free_x_num = 0;
            //free_x_num记录自由变元个数，如果超过1个，则解集无穷，无法求解.

            for(j=0 ; j<var ; j++){
                if(a[i][j] != 0 && free_x[j]) free_x_num++ , free_index = j;
            }

            if(free_x_num > 1) continue; //无法求解出确定变元

            tmp = a[i][var];
            for(j=0 ; j<var ; j++)
                if(a[i][j] != 0 && j != free_index)
                    tmp -= a[i][j]*x[j];
            x[free_index] = tmp / a[i][free_index];//求出该变元
            free_x[free_index] = false;//该变元确定
        }
        return var - k;//自由变元有var-k个
    }

    //唯一解的情况：在var*(var+1) 的增广矩阵中形成严格的上三角阵
    //计算x[n-1] , x[n-2] .... x[0] 回代计算，先得到后面的值，推到前面的值
    for(int i=var-1 ; i>=0 ; i--){
        tmp = a[i][var];
        for(int j=i+1 ; j<var ; j++){
            if(a[i][j] != 0 ) tmp -= a[i][j]*x[j];
        }
       // if(tmp % a[i][i] != 0) return -2; //说明浮点数存在，无整数解，但是有小数解
        x[i] = tmp / a[i][i];
    }
    return 0;
}