高斯消元法求解行列式

转自：http://blog.ac521.org/?p=53

利用高斯消元法形求解行列式的值，高斯消元就是这样的一个过程。
我们都知道行列式的代数性质（注意，第n+1列存储本方程等号右侧的系数）。这样，把第一行的要消去的元的系数和下面几行的相应的元的系数通过放缩化成一致（由于是实数，只放缩其中一个系数即可），进行减法即可消去该元。然后依次用第i行消去[i+1,n]的xi，直到只剩两项，即xn的系数和这时的等式右边的常数kn。只要用kn/xn就可以求出xn的近似解（在一定的数据范围内，double的精确度足够），然后依次回代，即可求出所有x值。
听着很玄乎，配合代码看就会很简单。
另外说一些关于精度的问题。实数类型不可避免的涉及精度的问题。显然，当n很大时，不断的进行除法运算会将近似更加近似，影响效果。这时候我们可以用整数进行运算，也就是找到两个数的最小公倍数进行减法就可以避免除法的近似。但是，整型会出现被0除的错误，需要特殊处理，并且当要消的元两个系数是两个大质数的话，行列式中的元素就会变得很大。如果再进行化简，就是又一重时间和空间的浪费，并且编程复杂度大大增高。所以，在一定的数据范围内，用实数完全可以实现高斯消元。
通过高斯消元把行列式转换成右上三角形的形势，再求解行列式（此时行列式的值就是对角线数值的相乘的积）

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
#include<iostream>
using namespace std;

const int n=4;
const double eps=1e-6;
double a[n][n]={
    3,  1, -1,  2,
    -5,  1,  3, -4,
    2,  0,  1, -1,
    1, -5,  3,  -3
}; //数组a是要求的行列式

bool zero(double a) //更精确的判断一个double类型的是否是0
{
    return a>-eps && a<eps;
}

double Gauss() {
    double mul,Result=1;
    int i,j,k,b[n];
    for(i=0;i<n;i++) b[i]=i;

    for(i=0;i<n;i++){
        if(zero(a[b[i]][i]))
            for(j=i+1;j<n;j++) //但是整行数据一个一个互换比较浪费，所以我在这里用b[i]数据当下表，然后更换一下b[i],b[j]所代表的值就达到的更换的效果
                if(!zero(a[b[j]][i])) {swap(b[i],b[j]); Result*=-1; break;  }
            Result*=a[b[i]][i];
            for(j=i+1;j<n;j++)
                if(!zero(a[b[j]][i])){
                    mul=a[b[j]][i]/a[b[i]][i];
                    for(k=i;k<n;k++)
                        a[b[j]][k]-=a[b[i]][k]*mul;
                }
    }
    return Result;
}

int main()
{
    cout<<Gauss()<<endl;
}