[数论] 高斯消元（整型和浮点型）

最新推荐文章于 2019-05-05 19:58:16 发布

同学少年

最新推荐文章于 2019-05-05 19:58:16 发布

阅读量3.5k

点赞数 1

分类专栏： ACM常用算法及模板整理文章标签：高斯消元

本文链接：https://blog.csdn.net/sr_19930829/article/details/41551567

版权

ACM常用算法及模板整理专栏收录该内容

20 篇文章 107 订阅

订阅专栏

高斯消元法：

数学上，高斯消元法（或译：高斯消去法）（英语：Gaussian Elimination），是线性代数中的一个算法，可用来为线性方程组求解，求出矩阵的秩，以及求出可逆方阵的逆矩阵。当用于一个矩阵时，高斯消元法会产生出一个“行梯阵式”。（来自维基百科）

构造如下方程：

a[0][0]*X0 + a[0][1] *X1 + a[0][2]*X2+...........a[0][n-1]*Xn-1 = a[0][n]
a[1][0]*X0 + a[1][1] *X1 + a[1][2]*X2+...........a[1][n-1]*Xn-1 ) = a[1][n]
..................................................
..................................................
a[m-1][0]*X0 + a[m-1][1] *X1 + a[m-1][2]*X2+...........a[m-1][n-1]*Xn-1 = a[m-1][n]

一共有m个方程，有n个未知量（X0，X1，...XN-1),未知量为所求，a[0....m-1][n]为常数。

在一些ACM题目中关键点就是构造方程，在求解概率期望的时候经常用到，找到题目中的状态递推方程。

常用E[k]表示从k节点到达目标节点还需要的期望的步数，那么目标节点p, E[p]=0，很关键。

构造出方程组，带入模板。

整型高斯消元模板：

//高斯消元模板  
#include <iostream>  
#include <stdio.h>  
#include <string.h>  
#include <string>  
#include <cmath>  
using namespace std;  
const int maxn=105;  
int equ,var; // 有equ个方程，var个变元。增广阵行数为equ, 分别为0到equ - 1，列数为var + 1，分别为0到var.  
int a[maxn][maxn];  
int x[maxn]; // 解集.  
bool free_x[maxn]; // 判断是否是不确定的变元.  
int free_num;  
void Debug(void)  
{  
    int i,j;  
    for(i=0;i<equ;i++)  
    {  
        for(j=0;j<var+1;j++)  
        {  
            cout<<a[i][j]<<" ";  
        }  
        cout<<endl;  
    }  
    cout<<endl;  
}  
inline int gcd(int a, int b)  
{  
    int t;  
    while (b!=0)  
    {  
        t=b;  
        b=a%b;  
        a=t;  
    }  
    return a;  
}  
inline int lcm(int a, int b)  
{  
    return a*b/gcd(a,b);  
}  
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)  
int Gauss(void)  
{  
    int i,j,k;  
    int max_r; // 当前这列绝对值最大的行.  
    int col; // 当前处理的列.  
    int ta,tb;  
    int LCM;  
    int temp;  
    int free_x_num;  
    int free_index;  
    // 转换为阶梯阵.  
    col=0; // 当前处理的列.  
    for(k=0;k<equ&&col<var;k++,col++)  
    { // 枚举当前处理的行.  
        // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)  
        max_r=k;  
        for(i=k+1;i<equ;i++)  
        {  
            if(abs(a[i][col])>abs(a[max_r][col]))  
                max_r=i;  
        }  
        if(max_r!=k)  
        { // 与第k行交换.  
            for(j=k;j<var+1;j++)  
                swap(a[k][j],a[max_r][j]);  
        }  
        if(a[k][col]==0)  
        { // 说明该col列第k行以下全是0了，则处理当前行的下一列.  
            k--; continue;  
        }  
        for(i=k+1;i<equ;i++)  
        { // 枚举要删去的行.  
            if (a[i][col]!=0)  
            {  
                LCM=lcm(abs(a[i][col]),abs(a[k][col]));  
                ta=LCM/abs(a[i][col]),tb=LCM/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;  
                }  
            }  
        }  
    }  
    //Debug();  
    // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).  
    for(i=k;i<equ;i++)  
    { // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.  
        if (a[i][col]!=0)  
            return -1;  
    }  
    // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.  
    // 且出现的行数即为自由变元的个数.  
    if(k<var)  
    {  
        // 首先，自由变元有var - k个，即不确定的变元至少有var - k个.  
        for (i=k-1;i>=0;i--)  
        {  
            // 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.  
            // 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.  
            free_x_num=0; // 用于判断该行中的不确定的变元的个数，如果超过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; // 无法求解出确定的变元.  
            // 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.  
            temp=a[i][var];  
            for(j=0;j<var;j++)  
            {  
                if(a[i][j]!=0&&j!=free_index)  
                temp-=a[i][j]*x[j];  
            }  
            x[free_index]=temp/a[i][free_index]; // 求出该变元.  
            free_x[free_index]=0; // 该变元是确定的.  
        }  
        return var-k; // 自由变元有var - k个.  
    }  
    // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.  
    // 计算出Xn-1, Xn-2 ... X0.  
    for (i=var-1;i>=0;i--)  
    {  
        temp=a[i][var];  
        for(j=i+1;j<var;j++)  
        {  
            if(a[i][j]!=0)   
                temp-=a[i][j]*x[j];  
        }  
        if(temp%a[i][i]!=0)   
            return -2; // 说明有浮点数解，但无整数解.  
        x[i]=temp/a[i][i];  
    }  
    return 0;  
}  
int main(void)  
{  
    int i, j;  
    while (scanf("%d %d",&equ,&var)!=EOF)  
    {  
        memset(a,0,sizeof(a));  
        memset(x,0,sizeof(x));  
        memset(free_x,1,sizeof(free_x)); // 一开始全是不确定的变元  
          
        for(i=0;i<equ;i++)//构造增广矩阵  
            for(j=0;j<var+1;j++)  
                scanf("%d",&a[i][j]);  
//        Debug();  
        free_num=Gauss();  
        if(free_num==-1) printf("无解!\n");  
        else if(free_num==-2) printf("有浮点数解，无整数解!\n");  
        else if(free_num>0)  
        {  
            printf("无穷多解! 自由变元个数为%d\n",free_num);  
            for(i=0;i<var;i++)  
            {  
                if(free_x[i]) printf("x%d 是不确定的\n",i+1);  
                else printf("x%d: %d\n",i+1,x[i]);  
            }  
        }  
        else  
        {  
            for(i=0;i<var;i++)  
                printf("x%d: %d\n",i+1,x[i]);  
        }  
        printf("\n");  
    }  
    return 0;  
}

浮点数高斯消元模板：

新的，用下面这个：

const int maxn=1002;
const double eps=1e-12;
double a[maxn][maxn];
int equ,var;//equ个方程,var个变量
double x[maxn];//解集
bool free_x[maxn];
int n;

int sgn(double x)
{
    return (x>eps)-(x<-eps);
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(0表示无解，1表示唯一解，大于1表示无穷解，并返回自由变元的个数)
int gauss()
{
    equ=n,var=n;//多少个方程，多少个变量
    int i,j,k;
    int max_r; // 当前这列绝对值最大的行.
    int col; // 当前处理的列.
    double temp;
    int free_x_num;
    int free_index;
    // 转换为阶梯阵.
    col=0; // 当前处理的列.
    memset(free_x,true,sizeof(free_x));
    for(k=0;k<equ&&col<var;k++,col++)
    {
        max_r=k;
        for(i=k+1;i<equ;i++)
        {
            if(sgn(fabs(a[i][col])-fabs(a[max_r][col]))>0)
                max_r=i;
        }
        if(max_r!=k)
        { // 与第k行交换.
            for(j=k;j<var+1;j++)
                swap(a[k][j],a[max_r][j]);
        }
        if(sgn(a[k][col])==0)
        { // 说明该col列第k行以下全是0了，则处理当前行的下一列.
            k--; continue;
        }
        for(i=k+1;i<equ;i++)
        { // 枚举要删去的行.
            if (sgn(a[i][col])!=0)
            {
                temp=a[i][col]/a[k][col];
                for(j=col;j<var+1;j++)
                {
                    a[i][j]=a[i][j]-a[k][j]*temp;
                }
            }
        }
    }

    for(i=k;i<equ;i++)
    {
        if (sgn(a[i][col])!=0)
            return 0;
    }
    if(k<var)
    {
        for(i=k-1;i>=0;i--)
        {
            free_x_num=0;
            for(j=0;j<var;j++)
            {
                if (sgn(a[i][j])!=0&&free_x[j])
                    free_x_num++,free_index=j;
            }
            if(free_x_num>1) continue;
            temp=a[i][var];
            for(j=0;j<var;j++)
            {
                if(sgn(a[i][j])!=0&&j!=free_index)
                    temp-=a[i][j]*x[j];
            }
            x[free_index]=temp/a[i][free_index];
            free_x[free_index]=0;
        }
        return var-k;
    }

    for (i=var-1;i>=0;i--)
    {
        temp=a[i][var];
        for(j=i+1;j<var;j++)
        {
            if(sgn(a[i][j])!=0)
                temp-=a[i][j]*x[j];
        }
        x[i]=temp/a[i][i];
    }
    return 1;
}

下面这个模板有些题目有问题......

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string.h>
#include <cmath>
#include <iomanip>
#include <algorithm>
using namespace std;

///浮点型高斯消元模板
const double eps=1e-12;
const int maxm=1000;///m个方程，n个变量
const int maxn=1000;
int m,n;
double a[maxm][maxn+1];///增广矩阵
bool free_x[maxn];///判断是否是不确定的变元
double x[maxn];///解集

int sign(double x)
{
    return (x>eps)-(x<-eps);
}
/**返回值：
-1 无解
0 有且仅有一个解
>=1 有多个解，根据free_x判断哪些是不确定的解
*/
int Gauss()
{
    int i,j;
    int row,col,max_r;
    m=n;///n个方程，n个变量的那种情况
    for(row=0,col=0;row<m&&col<n;row++,col++)
    {
        max_r=row;
        for(i=row+1;i<m;i++)///找到当前列所有行中的最大值(做除法时减小误差)
        {
            if(sign(fabs(a[i][col])-fabs(a[max_r][col]))>0)
                max_r=i;
        }
        if(max_r!=row)
        {
            for(j=row;j<n+1;j++)
                swap(a[max_r][j],a[row][j]);
        }
        if(sign(a[row][col])==0)///当前列row行以下全为0(包括row行)
        {
            row--;
            continue;
        }
        for(i=row+1;i<m;i++)
        {
            if(sign(a[i][col])==0)
                continue;
            double tmp=a[i][col]/a[row][col];
            for(j=col;j<n+1;j++)
                a[i][j]-=a[row][j]*tmp;
        }
    }
    for(i=row;i<m;i++)///col=n存在0...0,a的情况,无解
    {
        if(sign(a[i][col]))
            return -1;
    }
    if(row<n)///存在0...0,0的情况,有多个解,自由变元个数为n-row个
    {
        for(i=row-1;i>=0;i--)
        {
            int free_num=0;///自由变元的个数
            int free_index;///自由变元的序号
            for(j=0;j<n;j++)
            {
                if(sign(a[i][j])!=0&&free_x[j])
                    free_num++,free_index=j;
            }
            if(free_num>1)
                continue;///该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元
        ///只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的
            double tmp=a[i][n];
            for(j=0;j<n;j++)
            {
                if(sign(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 n-row;
    }
    ///有且仅有一个解,严格的上三角矩阵(n==m)
    for(i=n-1;i>=0;i--)
    {
        double tmp=a[i][n];
        for(j=i+1;j<n;j++)
        if(sign(a[i][j])!=0)
            tmp-=a[i][j]*x[j];
                x[i]=tmp/a[i][i];
    }
    return 0;
}///模板结束

int t,xx;

int main()
{
    cin>>t;
    while(t--)
    {
        cin>>n>>xx;
        memset(a,0,sizeof(a));
        for(int i=0;i<n;i++)
        {
            if(i==xx)
            {
                a[i][i]=1;
                a[i][n]=0;
                continue;
            }
            a[i][i]=1;
            a[i][n]=1;
            a[i][(i-1+n)%n]=-0.5;
            a[i][(i+1)%n]=-0.5;
        }
        Gauss();
        cout<<setiosflags(ios::fixed)<<setprecision(4)<<x[0]<<endl;
    }
    return 0;
}

整型高斯消元的一个题目：http://blog.csdn.net/sr_19930829/article/details/38959149

浮点高斯消元的一个题目：http://blog.csdn.net/sr_19930829/article/details/41551101