浮点数高斯消元:
bool Free[maxm];
int Gauss(int n, int m, double A[][maxm]){//A[i][0..m-1]系数矩阵 A[i][m]:常数矩阵
int row = 0;
for(int col = 0; row < n && col < m; row ++, col ++){
int mxr = row;
for(int i = row + 1; i < n; i ++){
if(fabs(A[i][col]) > fabs(A[mxr][col])) mxr = i;
}
if(mxr != row){
for(int j = 0; j <= m; j ++) swap(A[row][j], A[mxr][j]);
}
if(fabs(A[row][col]) < eps) { row --; continue; }
//将第col列出了第row行外全都消为0
for(int i = 0; i < n; i ++){
if(i == row) continue;
double xi = A[i][col] / A[row][col];
for(int j = m ;j >= col; j --){
A[i][j] -= xi * A[row][j];
}
}
}
//判断解的情况
//1.有唯一解
if(row == n && m == n){
for(int i = 0; i < n; i ++)
x[i] = A[i][m] / A[i][i];
return 0;
}
else{
//2.无解
for(int i = row; i < n; i ++)
if(fabs(A[i][m]) > eps){
return -1;
}
//3.多组解
memset(Free, true, sizeof(Free));
for(int i = row-1; i >= 0; i --){
int cnt = 0, id = -1;
//判断一行内自由变元的个数
for(int j = 0; j < m; j ++){
if(fabs(A[i][j]) > eps && Free[j]) { cnt++; id = j; } //该eps不可用0替代,否则会出问题
}
if(!cnt || cnt > 1) continue;
//一行只有一个自由变元,可解
double tmp = A[i][m];
for(int j = 0; j < m; j ++){
if(j != id && fabs(A[i][j]) > eps) tmp -= A[i][j] * x[j];
}
x[id] = tmp / A[i][id];
Free[id] = 0;
}
return m - row;//返回自由变元的个数
}
}