hihocoder #1195 高斯消元一

hihocoder对算法解释得很详细，就复制粘贴来了

首先我们要计算出上三角矩阵，也就是将方程组变为：

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] = d
                      e * x[3] = f
                             0 = 0

当我们在第一个式子中消去x[3]后，有a * x[1] + b * x[2] = g，显然x[1]和x[2]有无穷多种可能的取值。

然后贴一发模板

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iostream>
using namespace std;
#define db double
const int maxn=1008;
db a[maxn][maxn];
db x[maxn];
/*void Debug(int equ,int var)
{
    int i, j;
    for (i = 1; i <=equ; i++)
    {
        for (j = 1; j <= var + 1; j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}*/
int guess(int row,int col)
{
    int i,j,k;
    for(i=1;i<col;i++)
    {
        int r=i;
        for(j=row;j>i;j--)
        {
            if(fabs(a[j][i])>fabs(a[r][i]))
                r=j;
        }
        if(r==i&&fabs(a[i][i])<1e-7)    return -1;
        if(r!=i)    swap(a[r],a[i]);
        for(j=i+1;j<=row;j++)
        {
            for(k=col;k>i;k--)
                a[j][k]-=a[j][i]/a[i][i]*a[i][k];
            a[j][i]=0;
        }
    }
//    Debug(row,col);
    for(i=col-1;i<=row;i++)
    {
        for(j=1;j<col;j++)
            if(fabs(a[i][j])>1e-6)
                break;
        if(j==col&&fabs(a[i][col])>1e-6)
            return 0;
    }
    for(i=col-1;i>0;i--)
    {
        for(j=i+1;j<col;j++)
            a[i][col]-=a[i][j]*x[j];
        x[i]=a[i][col]/a[i][i];
    }
    return 1;
}
int main()
{
    int n,m,i,j;
    scanf("%d%d",&n,&m);
    for(i=1;i<=m;i++)
        for(j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]); 
    int ans=guess(m,n+1);
    if(ans==-1)    printf("Many solutions
");
    else if(ans==0)    printf("No solutions
");
    else    
    {
        for(i=1;i<=n;i++)
            printf("%d
",(int)(x[i]+0.5));
    }
    return 0;
}