[模板]高斯消元

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define WR WinterRain
using namespace std;
const int WR=1010;
const double eps=1e-6;//焯!!!!!!!!精度不是int!!!!!!! 
int equ,var;//有equ个方程,var个变元
double a[WR][WR];//增广矩阵 
double x[WR];//存解集 
bool fre[WR];//标记是否为自由元,即可随意取值 
bool f=true;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<3)+(s<<1)+ch-48;
        ch=getchar();
    }
    return s*w;
}
void gauss(){
    /*转化成阶梯矩阵*/
    int row=1;
    for(int col=1;col<=var;col++){
        //row:当前处理的行
        //col:当前处理的列 
        int mxrow=row;//寻找当前列系数绝对值最大的行
        for(int i=row+1;i<=equ;i++){//没有必要去上面找,上面的都已经处理好了 
            if(fabs(a[i][col])>fabs(a[mxrow][col])){
                mxrow=i;
            }
        }
        if(fabs(a[mxrow][col])<eps){//col列从第row行以下全都为0,处理当前行的下一列 
        //    printf("111");
            continue;
        }
        //交换 
        for(int i=row;i<=var+1;i++){//手动模拟高斯消元可以发现
            //前面的几项(从1到row-1)都为0了,不用交换 
            swap(a[row][i],a[mxrow][i]);
        }
        for(int i=row+1;i<=equ;i++){//手模一遍相当于消去该行的第col列的元 
            if(fabs(a[i][col])>eps){
                double tmp=a[i][col]/a[row][col];//乘数 
                for(int j=col+1;j<=var+1;j++){
                    a[i][j]-=a[row][j]*tmp;
                }
                //a[i][col]=0.0;
            }
        }
        row++;
    }
    //printf("%d/n",row);
    //for(int i=1;i<=equ;i++){
    //    for(int j=1;j<=equ+1;j++){
    //        printf("%.2lf ",a[i][j]);
    //    }
    //    printf("/n");
    //}
    if(row<=equ){
        while(row<=equ){
            if(a[row++][var+1]!=0){
                printf("-1");
                f=false;
                return;
            }
        }
        printf("0");
        f=false;
        return;
    }
    for(int i=equ;i>=1;i--){
        double tmp=a[i][var+1];
        for(int j=i+1;j<=var;j++){
            tmp-=a[i][j]*x[j];
        }
        x[i]=tmp/a[i][i];
        //printf("%d %lf/n",i,x[i]);
    } 
}
int main(){
    equ=var=read();//这里让var=equ(看题面) 
    for(int i=1;i<=equ;i++){
        for(int j=1;j<=equ+1;j++){
            int t=read();
            a[i][j]=(double)t;
        }
    }
    gauss();
    if(!f) return 0;
    for(int i=1;i<=var;i++){
        if(fabs(x[i])<eps) printf("x%d=0/n",i);
        else printf("x%d=%.2lf/n",i,x[i]);
    }
    return 0;
}