本文主要是介绍hiho一下第56周 高斯消元，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
小Ho：<吧唧><吧唧><吧唧>
小Hi：小Ho，你还吃呢。想好了么?
小Ho：肿抢着呢(正想着呢)...<吞咽>...我记得这个问题上课有提到过，应该是一元一次方程组吧。
我们把每一件商品的价格看作是x[1]..x[n]，第i个组合中第j件商品数量记为a[i][j]，其价格记作y[i]，则可以列出方程式：
a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1] a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]... a[m][1] * x[1] + a[m][2] * x[2] + ... + a[m][n] * x[n] = y[m]
我们可以对方程组进行3种操作而不改变方程组的解集：
1. 交换两行。
2. 把第i行乘以一个非0系数k。即对于j = 1..n， 令a[i][j] = k*a[i][j], y[i]=k*y[i]
3. 把第p行乘以一个非0系数k之后加在第i行上。即对于j=1..n， 令a[i][j] = a[i][j]+k*a[p][j]，y[i]=y[i]+k*p[i]
以上三个操作叫做初等行变换。
我们可以使用它们，对这个方程组中的a[i][j]进行加减乘除变换，举个例子：
a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1] 式子(1) a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2] 式子(2)
我们可以通过 式子(1) - 式子(2) * (a[1][1] / a[2][1])，将第1行第1列的a[1][1]变换为0。
对整个方程组进行多次变换之后，可以使得a[i][j]满足：
a[i][j] = 1 (i == j) a[i][j] = 0 (i != j)
则整个方程组变成了：
x[1] = y'[1] x[2] = y'[2] ... x[n] = y'[n] 0 = y'[n + 1] 0 = y'[n + 2] ... 0 = y'[m]
这样的话，y'[1] .. y'[n]就是我们要求的x[1]..x[n]
小Hi：挺不错的嘛，继续？
小Ho：好，关于如何变换，我们可以利用一个叫高斯消元的算法。高斯消元分成了2个步骤：
首先我们要计算出上三角矩阵，也就是将方程组变为：
a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y'[1]0 * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y'[2]0 * x[1] + 0 * x[2] + ... + a[3][n] * x[n] = y'[3]...0 * x[1] + 0 * x[2] + ... + a[n][n] * x[n] = y'[n]0 * x[1] + 0 * x[2] + ... + 0 * x[n] = y'[n + 1]...0 * x[1] + 0 * x[2] + ... + 0 * x[n] = y'[m]
也就是通过变换，将所有a[i][j](i>j)变换为0。同时要保证对角线上的元素a[i][i]不为0。
方法也很见简单，从第1行开始，我们利用当前行第i列不为0，就可以通过变换将i+1..M行第一列全部变换为0，接着对于第2行，我们用同样的方法将第3..M行第2列也变换为0...不断重复直到第n行为止。
假如计算到第i行时，第i列已经为0，则我们需要在第i+1..M行中找到一行第i列不为0的行k，并交换第i行和第k行，来保证a[i][i] != 0。但这时候还有可能出现一个情况，就是第i..M行中的i列均为0，此时可以判定，该方程组有多解。
当得到上三角矩阵后，就可以从第n行开始逆推，一步一步将a[i][j](i<j)也变换为0.
因为第n行为a[n][n] * x[n] = y'[n]，则x[n] = y'[n] / a[n][n]。
第n-1行为a[n-1][n-1] * x[n - 1] + a[n][n] * x[n] = y'[n - 1]。我们将得到的x[n]代入，即可计算出x[n-1]。
同样的依次类推就可以得到所有的x[1]..x[n]。
而对于多解和无解的判定：
当在求出的上三角矩阵中出现了 a[i][1] = a[i][2] = ... = a[i][n] = 0, 但是y'[i] != 0时，产生了矛盾，即出现了无解的情况。
而多解的证明如下：
假设n=3,m=3，而我们计算出了上三角矩阵为：
a * x[1] + b * x[2] + c * x[3] = de * x[3] = f0 = 0
当我们在第一个式子中消去x[3]后，有a * x[1] + b * x[2] = g，显然x[1]和x[2]有无穷多种可能的取值。
小Hi：既然小Ho你都已经把整个算法讲了，那么我就只能给出伪代码了：
// 处理出上三角矩阵 For i = 1..NFlag ← FalseFor j = i..M // 从第i行开始，找到第i列不等于0的行jIf a[j][i] != 0 ThenSwap(j, i) // 交换第i行和第j行Flag ← TrueBreakEnd IfEnd For// 若无法找到，则存在多个解If (not Flag) ThenmanySolutionsFlag ← True // 出现了秩小于N的情况continue;End If// 消除第i+1行到第M行的第i列For j = i+1 .. Ma[j][] ← a[j][] - a[i][] * (a[j][i] / a[i][i])b[j] ← b[j] - b[i] * (a[j][i] / a[i][i])End For End For // 检查是否无解，即存在 0 = x 的情况 For i = 1..MIf (第i行系数均为0 and (b[i] != 0)) Thenreturn "No solutions"End If End For// 判定多解 If (manySolutionsFlag) Thenreturn "Many solutions" End If// 此时存在唯一解 // 由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数 // 解析来从第N行开始处理每一行的解 For i = N .. 1// 利用已经计算出的结果，将第i行中第i+1列至第N列的系数消除For j = i + 1 .. Nb[i] ← b[i] - a[i][j] * value[j]a[i][j] ← 0End Forvalue[i] ← b[i] / a[i][i] End For
那最后能够拜托你实现一下这个算法么？
小Ho：没问题，等我吃完这包薯片就去！
/************************************************ Author: fisty* Created Time: 2015-08-19 13:22:23* File Name : 56_Gauss.cpp*********************************************** */
#include <iostream>
#include <cstring>
#include <deque>
#include <cmath>
#include <queue>
#include <stack>
#include <list>
#include <map>
#include <set>
#include <string>
#include <vector>
#include <cstdio>
#include <bitset>
#include <algorithm>
using namespace std;
#define Debug(x) cout << #x << " " << x <<endl
#define Memset(x, a) memset(x, a, sizeof(x))
const int INF = 0x3f3f3f3f;
typedef long long LL;
typedef pair<int, int> P;
#define FOR(i, a, b) for(int i = a;i < b; i++)
#define lson l, m, k<<1
#define rson m+1, r, k<<1|1const double eps = 1e-8;
const int Max_M = 1000; ///m个方程,n个变量
const int Max_N = 500;
int m, n;
double Aug[Max_M][Max_N+1]; ///增广矩阵
bool free_x[Max_N]; ///判断是否是不确定的变元
double x[Max_N]; ///解集 int sign(double x){ return (x>eps) - (x<-eps);} /** 返回值: -1 无解 0 有且仅有一个解 >=1 有多个解,根据free_x判断哪些是不确定的解 */
int Gauss()
{ int i,j; int row,col,max_r; for(row=0,col=0; row<m&&col<n; row++,col++) { max_r = row; for(i = row+1; i < m; i++) ///找到当前列所有行中的最大值(做除法时减小误差) { if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0) max_r = i; } if(max_r != row) ///将该行与当前行交换 { for(j = row; j < n+1; j++) swap(Aug[max_r][j],Aug[row][j]); } if(sign(Aug[row][col])==0) ///当前列row行以下全为0(包括row行) { row--; continue; } for(i = row+1; i < m; i++) { if(sign(Aug[i][col])==0) continue; double ta = Aug[i][col]/Aug[row][col]; for(j = col; j < n+1; j++) Aug[i][j] -= Aug[row][j]*ta; } } for(i = row; i < m; i++) ///col=n存在0...0,a的情况,无解 { if(sign(Aug[i][col])) return -1; } if(row < n) ///存在0...0,0的情况,有多个解,自由变元个数为n-row个 { for(i = row-1; i >=0; i--) { int free_num = 0; ///自由变元的个数 int free_index; ///自由变元的序号 for(j = 0; j < n; j++) { if(sign(Aug[i][j])!=0 && free_x[j]) free_num++,free_index=j; } if(free_num > 1) continue; ///该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元 ///只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的 double tmp = Aug[i][n]; for(j = 0; j < n; j++) { if(sign(Aug[i][j])!=0 && j!=free_index) tmp -= Aug[i][j]*x[j]; } x[free_index] = tmp/Aug[i][free_index]; free_x[free_index] = false; } return n-row; } ///有且仅有一个解,严格的上三角矩阵(n==m) for(i = n-1; i >= 0; i--) { double tmp = Aug[i][n]; for(j = i+1; j < n; j++) if(sign(Aug[i][j])!=0) tmp -= Aug[i][j]*x[j]; x[i] = tmp/Aug[i][i]; } return 0;
} void init(){memset(Aug,0.0,sizeof(Aug)); memset(x,0.0,sizeof(x)); memset(free_x,1,sizeof(free_x)); ///都是不确定的变元
}
int main(){scanf("%d%d", &n, &m);init();int ans;for(int i = 0;i < m; i++){for(int j = 0;j <= n; j++){scanf("%lf", &Aug[i][j]);}}ans = Gauss();if(ans == -1){puts("No solutions");}else if(ans >= 1){puts("Many solutions");}else{for(int i = 0;i < n; i++){printf("%d\n", int(x[i] + eps));}}
}
这篇关于hiho一下第56周 高斯消元的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
