高斯消元——浮点数模板

#include <iostream>
#include <stdio.h>
#include <string.h>
#define maxn 210
#define CL(a,num) memset((a),(num),sizeof(a))
#define iabs(x)  ((x) > 0 ? (x) : -(x))
#define EPS 1e-8
using namespace std;
int equ,var;//方程个数和自由元的个数
double a[maxn][maxn];
double x[maxn],free_x[maxn];//解集

void Debug()
{
    int i,j;
    for(i=0;i<equ;i++){
        for(j=0;j<var+1;j++)
            printf("%.2lf ",a[i][j]);
        cout<<endl;
    }
}

int zero(double x)
{
    if( iabs(x)<EPS ) return 1;
    return 0;
}
int Guass(){
    int i,j,k,col;
    CL(x,0);CL(free_x,1);//清空解集
    
    for(k=0,col=0;k<equ && col<var ;k++,col++){//枚举行和列
        //    cout<<k<<" "<<col<<endl;
        int max_r = k;
        for(i=k+1;i<equ;i++)
            if(iabs(a[i][col]) > iabs(a[max_r][col])) max_r = i;
        if(max_r != k)//交换
            for(i=k;i<var+1;i++)
                swap(a[k][i],a[max_r][i]);
        
        if( zero(a[k][col]) ) {
            k--; //为什么有k--;模拟下这样可以减的多点
            continue;
        }
        
        for(i=k+1;i<equ;i++){
            if(!zero(a[i][col])){
                //    int lcm = LCM(a[k][col],a[i][col]);
                //    int ta = lcm/iabs(a[i][col]);
                //    int tb = lcm/iabs(a[k][col]);
                //    如果读入的矩阵是整数可按照上面的写，否则得按照下面的写
                double tb = iabs(a[i][col])/iabs(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] - tb*a[k][j];
            }
        }
    }
    //    Debug();
    //1.无解的情况出现(0,0,0,0,……a)这样的行且a!=0
    for(i=k;i<equ;i++){
        if(!zero(a[i][col]))
        {
            return -1;
        }
    }
    //2.无穷解的情况，再var*(var+1)的增广矩阵中出现(0,0,0……0)这样的行
    //出现的行数便是自由变元的个数
    if(k<var){
        int num = 0,freeidx=0;
        for(i=k-1;i>=0;i--){
            num=0;
            //用于判断该行的不确定的变元的个数，如果超过1则无法求解
            //它们仍然为不确定的变元
            double tmp = a[i][var];
            for(j=0;j<var;j++)
                if(!zero(a[i][j]) && free_x[j]){
                    num++;
                    freeidx = j;
                }
            if(num>1) continue;
            tmp=a[i][var];
            for(j=0;j<var;j++){
                if(!zero(a[i][j]) && j != freeidx) tmp -= a[i][j]*x[j];
            }
            x[freeidx] = tmp/a[i][freeidx];
            free_x[freeidx]=0;
        }
        
        return var-k;
    }
    //3.唯一解的情况
    for(i=k-1;i>=0;i--)
    {
        double tmp=a[i][var]*1.0;
        for(j=i+1;j<var;j++)
            tmp-=a[i][j]*x[j];
        x[i] = tmp/(a[i][i]*1.0);
    }
    return 0;
}

//int main()
//{
//    int t,co=0;
//    scanf("%d",&t);
//    while(t--)
//    {
//        int i,j,k;
//        scanf("%d",&equ);
//        var=equ;
//        for(i=0;i<equ;i++)
//            for(j=0;j<var+1;j++){
//                scanf("%lf",&a[i][j]);
//            }
//        //Debug();
//        //cout<<endl;
//        int free_num=Guass();
//        co++;
//        printf("Case #%d:",co);
//        if(free_num==0)
//        {
//            printf("%.2lf",(iabs(x[0])<EPS?EPS:x[0]));
//            for(i=1;i<var;i++) printf(" %.2lf",(iabs(x[i])<EPS?EPS:x[i]));
//            cout<<endl;
//        }else cout<<"Can't solve it."<<endl;
//        //无求无穷多解的情况
//    }
//    return 0;
//}