#1195 : 高斯消元·一
-
2 2 2 1 5 1 2 4
样例输出
-
2 1
描述
小Ho:喂不得了啦,那边便利店的薯片半价了!
小Hi:啥?!
小Ho:那边的便利店在打折促销啊。
小Hi:走走走,赶紧去看看=v=
于是小Hi和小Ho来到了便利店。
老板为了促销,推出了组合包的形式,将不同数量的各类商品打包成一个组合,顾客可以选择组合进行购买。比如2袋薯片,1听可乐的组合只要5元,而1袋薯片,2听可乐的组合只要4元。
通过询问老板,小Hi和小Ho知道:一共有N种不同的商品和M种不同的商品组合;每一个组合的价格等于组合内商品售价之和,一个组合内同一件商品不会超过10件。
小Hi:这样算下来的话,一听可乐就是1元,而一包薯片是2元。小Ho,如果你知道所有的组合情况,你能分别算出每一件商品单独的价格么?
小Ho:当然可以了,这样的小问题怎么能难到我呢?
输入
第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个商品单独的售价。数据保证如果存在唯一解,那么解一定恰好是非负整数解。
提示:高斯消元
小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 manySolutionsFlag ← True // 出现了秩小于N的情况 continue; 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 // 判定多解 If (manySolutionsFlag) Then return "Many solutions" End If // 此时存在唯一解 // 由于每一行都比前一行少一个系数,所以在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:没问题,等我吃完这包薯片就去!
呜呜,我开始一直用的longlong呀-----WA了一页了-----
用longlong相消时,最小公倍数可能超long long ----
附下自认为完整的long long版---
#include<cstdio>
#include<cstring>
#include<stack>
#include<cmath>
#include<queue>
#include<algorithm>
using namespace std;
#define LL long long
LL shu[1200][600],ans[1200];
LL gcd(LL a,LL b)
{
if (a<b)
{
LL c=a;
a=b;
b=c;
}
if (b==0) return a;
return gcd(b,a%b);
}
LL gg(LL a,LL b)
{
int aa,bb;
aa=bb=0;
if (a<0) a=-a,aa=1;
if (b<0) b=-b,bb=1;
LL yu=gcd(a,b);
a/=yu;
if (aa) a=-a;
if (bb) a=-a;
return a;
}
int main()
{
int n,m,i,j;//printf("%lld %lld %lld %lld\n",gcd(0,2),gcd(2,0),gcd(1,2),gcd(0,0));
LL yu;bool fafe[2];
fafe[1]=fafe[0]=false;
scanf("%d%d",&n,&m);
for (i=1;i<=m;i++)
for (j=1;j<=n+1;j++)
scanf("%lld",&shu[i][j]);
//输入完成===
stack <int > que;
for (i=1;i<=m;i++)//化简---消零---
{
yu=abs(shu[i][1]);
for (j=2;j<=n+1;j++)
yu=gcd(yu,abs(shu[i][j]));
if (yu)
{
if (yu!=1)
for (j=1;j<=n+1;j++)
shu[i][j]/=yu;
}
else
que.push(i);
}
while (!que.empty())
{
int k=que.top();
que.pop();
for (int pp=1;pp<=n+1;pp++)
shu[k][pp]=shu[m][pp];
m--;
}
int lie=1,op,c;
for (i=1;i<=n+1&&lie<m;i++)
{
op=0;//找一个近0首项
for (j=lie;j<=m;j++)
if (shu[j][i])
{
if (!op)
op=j;
else if (abs(shu[j][i])<abs(shu[op][i]))
op=j;
}
if (!op)//这一列全0--
continue;
if (op!=lie)//交换lie行和op行
{
for (int k=1;k<=n+1;k++)
{
c=shu[op][k];
shu[op][k]=shu[lie][k];
shu[lie][k]=c;
}
}
for (int k=lie+1;k<=m;k++)//用lie行来消下面的行
{
if (!shu[k][i]) continue;
yu=gg(shu[lie][i],shu[k][i]);
LL kop=shu[k][i]*yu/shu[lie][i];//通分
for (int kl=i;kl<=n+1;kl++)//消--
shu[k][kl]=shu[k][kl]*yu-kop*shu[lie][kl];
yu=abs(shu[k][i]);//化简
for (int kl=i+1;kl<=n+1;kl++)
yu=gcd(yu,abs(shu[k][kl]));
if (yu)
{
if (yu!=1)
{
for (int kl=i;kl<=n+1;kl++)
shu[k][kl]/=yu;
}
}
else
que.push(k);
}
lie++;
while (!que.empty())//消空行
{
int k=que.top();
que.pop();
for (int pp=1;pp<=n+1;pp++)
shu[k][pp]=shu[m][pp];
m--;
}
}
// for (i=1;i<=m;i++)
// {
// for (j=1;j<=n+1;j++)
// printf("%lld ",shu[i][j]);
// printf("\n");
// }
fafe[1]=true;//0=C--------无解的情况
for (i=1;i<=n;i++)
if (shu[m][i])
fafe[1]=false;
if (fafe[1])
printf("No solutions\n");
else if (m!=n)
printf("Many solutions\n");
else
{
for (j=n;j>0;j--)
{
ans[j]=shu[j][n+1];
for (int k=n;k>j;k--)
ans[j]-=ans[k]*shu[j][k];
ans[j]/=shu[j][j];
}
for (i=1;i<=n;i++)
printf("%lld\n",ans[i]);
}
return 0;
}
用double的需要自设一个近0值---10^-n 次方---
代码:
#include<cstdio>
#include<cstring>
#include<stack>
#include<cmath>
#include<queue>
#include<algorithm>
using namespace std;
#define LL long long
double shu[1200][600],ans[1200];
int cinn[1200];
int main()
{
int n,m,i,j;
LL yu;bool fafe[2];
while (~scanf("%d%d",&n,&m))
{
fafe[1]=fafe[0]=false;
double ling=0.000001;
for (i=1;i<=m;i++)
for (j=1;j<=n+1;j++)
scanf("%lf",&shu[i][j]);
//输入完成===
stack <int > que;
bool yy,xx;
double xie;
int lie=1,op;//lie表示我们要用第lie行去消下面的行
double c;
for (i=1;i<=n;i++)
{
op=0;//找一个近0首项
for (j=lie;j<=m;j++)
if (fabs(shu[j][i])>ling)
{
op=j;
break;
}
if (!op)//这一列全0或行数不够--不是无解就是多解-.-然后i后移---用lie行消i+1列
{
fafe[0]=true;
continue;
}
if (op!=lie)//交换lie行和op行
{
for (int k=1;k<=n+1;k++)
{
c=shu[op][k];
shu[op][k]=shu[lie][k];
shu[lie][k]=c;
}
}
for (int k=lie+1;k<=m;k++)//用lie行来消下面的行
{
xie=shu[k][i]/shu[lie][i];//这里精良不要用b=b*k-a;---会使后项越变越大-.-
for (int mm=i;mm<=n+1;mm++)
shu[k][mm]-=shu[lie][mm]*xie;
}
for (int k=m;k>lie;k--)//消0行-.-这一步也可以无视
{
xx=true;
for (j=1;j<=n+1;j++)
if (fabs(shu[k][j])>ling)
{
xx=false;
break;
}
if (xx)
{
if (k!=m)
for (j=1;j<=n+1;j++)
shu[k][j]=shu[m][j];
m--;
}
}
lie++;
}
for (i=1;i<=m;i++)//无解的情况--
{
if (fabs(shu[i][n+1])>ling)
fafe[1]=true;
for (j=1;j<=n;j++)
if (fabs(shu[i][j])>ling)
{
fafe[1]=false;break;
}
if (fafe[1]) break;
}
if (fafe[1])
{
printf("No solutions\n");
continue;
}
else if (fafe[0])
printf("Many solutions\n");
else
{
for (j=n;j>0;j--)
{
ans[j]=shu[j][n+1];
for (int k=n;k>j;k--)
ans[j]-=ans[k]*shu[j][k];
ans[j]/=shu[j][j];
}
for (i=1;i<=n;i++)
cinn[i]=(int)(ans[i]+0.5);
for (int i=1;i<=n;i++)
printf("%d\n",cinn[i]);
}
}
return 0;
}