hiho一下 第五十六周



题目1 : 高斯消元·一

时间限制: 10000ms

单点时限: 1000ms

内存限制: 256MB

描述

小Ho：喂不得了啦，那边便利店的薯片半价了!

小Hi：啥?!

小Ho：那边的便利店在打折促销啊。

小Hi：走走走，赶紧去看看=v=

于是小Hi和小Ho来到了便利店。

老板为了促销，推出了组合包的形式，将不同数量的各类商品打包成一个组合，顾客可以选择组合进行购买。比如2袋薯片，1听可乐的组合只要5元，而1袋薯片，2听可乐的组合只要4元。

通过询问老板，小Hi和小Ho知道：一共有N种不同的商品和M种不同的商品组合；每一个组合的价格等于组合内商品售价之和，一个组合内同一件商品不会超过10件。

小Hi：这样算下来的话，一听可乐就是1元，而一包薯片是2元。小Ho，如果你知道所有的组合情况，你能分别算出每一件商品单独的价格么？

小Ho：当然可以了，这样的小问题怎么能难到我呢？

   

提示：高斯消元

 

          ×          

提示：高斯消元

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

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

小Hi：既然小Ho你都已经把整个算法讲了，那么我就只能给出伪代码了：

// 处理出上三角矩阵
For i = 1..N
    Flag ← False
    For j = i..M                // 从第i行开始，找到第i列不等于0的行j
        If a[j][i] != 0 Then
            Swap(j, i)          // 交换第i行和第j行
            Flag ← True
            Break
        End If
    End For
    // 若无法找到，则存在多个解
    If (not Flag) Then
        return "Many solutions"
    End If
    // 消除第i+1行到第M行的第i列
    For j = i+1 .. M
        a[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..M
    If (第i行系数均为0 and (b[i] != 0)) Then
        return "No solutions"
    End If
End For

// 此时存在唯一解
// 由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数
// 解析来从第N行开始处理每一行的解
For i = N .. 1
    // 利用已经计算出的结果，将第i行中第i+1列至第N列的系数消除
    For j = i + 1 .. N
        b[i] ← b[i] - a[i][j] * value[j]
        a[i][j] ← 0
    End For
    value[i] ← b[i] / a[i][i]
End For
		

那最后能够拜托你实现一下这个算法么？

小Ho：没问题，等我吃完这包薯片就去！

          Close        

输入

第1行：2个正整数，N,M。表示商品的数量N，组合的数量M。1≤N≤500, N≤M≤2*N

第2..M+1行：N+1个非负整数，第i+1行第j列表示在第i个组合中，商品j的数量a[i][j]。第i+1行第N+1个数表示该组合的售价c[i]。0≤a[i][j]≤10, 0≤c[i]≤10^9

输出

若没有办法计算出每个商品单独的价格，输出"No solutions"

若可能存在多个不同的结果，输出"Many solutions"

若存在唯一可能的结果，输出N行，每行一个非负整数，第i行表示第i个商品单独的售价。数据保证如果存在唯一解，那么解一定恰好是非负整数解。

样例输入

2 2
2 1 5
1 2 4

样例输出

2
1

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
double a[1050][550],b[1050],ans[550];//a为每个商品的数量，b为总价，ans为最后的单价  
const double fbs=1e-6;
int n,m;
int solve()
{
    for(int i=1;i<=n;i++)
    {
        int flag=0;
        for(int j=i;j<=m;j++) //从第i行开始，找到第i列不等于0的行j 
        {
            if(a[j][i]!=0)
            {
                for(int k=1;k<=n;k++)
                    swap(a[i][k],a[j][k]); //交换第i行和第j行 
                swap(b[i],b[j]);
                flag=1;
            }
            if(flag) break;
        }
        if(!flag)   //若无法找到，则存在多个解  
        {
            printf("Many solutions\n");
            return 0;
        }
        //消除第i+1行到第M 行的第i列  
        for(int j=i+1;j<=m;j++)
        {
            double div=a[j][i]/a[i][i];
            for(int k=i;k<=n;k++)
            {
                a[j][k]=a[j][k]-a[i][k]*div;
            }
            b[j]=b[j]-b[i]*div;
        }
    }
    //检查是否无解，存在0=x的情况 
    for(int i=1;i<=m;i++)   //注意i<=m,之前写成i<=n
    {
        int flag=0;
        for(int j=1;j<=n;j++)
        {
            if(a[i][j]<-fbs || a[i][j]>fbs)
            {
                 flag=1;
                 break;
            }
        }
        if(flag==0 && (b[i]>fbs || b[i]<-fbs))
        {
           printf("No solutions\n");
           return 0;
        }
    }
    //存在唯一解，  
    //由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数  
    //解析从第Ｎ航开始处理每一行的解  
    for(int i=n;i>=1;i--)
    {
        for(int j=i+1;j<=n;j++)
        {
            b[i]=b[i]-a[i][j]*ans[j];
            a[i][j]=0;
        }
        ans[i]=b[i]/a[i][i];
    }
    return 1;
}


int main()
{
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        for(int i=1;i<=m;i++)
        {
            for(int j=1;j<=n;j++)
                scanf("%lf",&a[i][j]);
            scanf("%lf",&b[i]);
        }
        if(solve())
        {
            for(int i=1;i<=n;i++)
                printf("%d\n",(int) (ans[i]+0.5));
        }

    }
    return 0;
}