高斯_约当消元(浮点)

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
#define MAXN 400
const double eps=1e-8;

int n,m;
double a[MAXN+10][MAXN+10],x[MAXN+10];
bool free_x[MAXN+10];

void Debug();

void read()
{
    int x;
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m+1;j++){
            scanf("%d",&x);
            a[i][j]=1.0*x;
        }
    }
}
void GJ_elimination(int equ,int var,int &row,int &col)
{
    int mx;
    for(row=col=1;row<=equ&&col<=var;row++,col++){
        mx=row;
        for(int i=row+1;i<=equ;i++){
            if(fabs(a[i][col])>fabs(a[mx][col]))
                mx=i;
        }
        if(mx!=row){
            for(int j=1;j<=var+1;j++)
                swap(a[mx][j],a[row][j]);
        }
        if(fabs(a[row][col])<eps){
            row--;
            continue;
        }
        for(int i=1;i<=equ;i++){
            if(i==row||fabs(a[i][col])<eps)
                continue;
            for(int j=var+1;j>=col;j--)
                a[i][j]-=a[i][col]/a[row][col]*a[row][j];
        }
    }
}
int GJ_Judge(int equ,int var,int row,int col)
{
    //-1:no solution
    for(int i=row;i<=equ;i++)
        if(fabs(a[i][var+1])>eps)
            return -1;
    //>0:不定解
    if(row<=var){
        memset(free_x,0,sizeof free_x);
        for(int i=row-1;i>=1;i--){
            int freex=0,tar;
            for(int j=1;j<=var;j++){
                if(fabs(a[i][j])>eps&&!free_x[j])
                    freex++,tar=j;
            }
            if(freex>1) continue;
            double tmp=0.0;
            for(int j=1;j<=var;j++)
                if(fabs(a[i][j])>eps&&j!=tar)
                    tmp+=a[i][j]*x[j];
            x[tar]=(a[i][var+1]-tmp)/a[i][tar];
            free_x[tar]=true;
        }
        return var-row+1;
    }
    //0:唯一解
    for(int i=row-1;i>=1;i--){ // Attention!
        double tmp=0.0;
        for(int j=i+1;j<=var;j++)
            if(fabs(a[i][j])>eps)
                tmp+=a[i][j]*x[j];
        x[i]=(a[i][var+1]-tmp)/a[i][i];
    }
    return 0;
}
void Debug()
{
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m+1;j++)
            printf("%.0lf ",a[i][j]);
            puts("");
        }
        puts("");
}
void Gauss_Jordan(int equ,int var)
{
    int row,col;
    memset(x,0,sizeof x);
    GJ_elimination(equ,var,row,col);
    int p=GJ_Judge(equ,var,row,col);
    if(p==-1)
        printf("No solution\n");
    else if(p==0){
        for(int j=1;j<=var;j++)
            printf("X[%d] = %.4lf\n",j,x[j]+eps);
    }
    else{
        printf("Multiple solution, free_num: %d\n",p);
        for(int j=1;j<=var;j++){
            if(!free_x[j])
                printf("X[%d] not determined\n",j);
            else
                printf("X[%d] = %.4lf\n",j,x[j]+eps);
        }
    }
}
int main()
{
    read();
    Gauss_Jordan(n,m);
}