hiho一下 第五十六周 #1195 : 高斯消元·一 【高斯消元--】

最新推荐文章于 2019-07-18 08:05:10 发布
最新推荐文章于 2019-07-18 08:05:10 发布
阅读量510 收藏
点赞数
分类专栏： hihocoder 数学
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/leibniz_zhang/article/details/53283923
版权

#1195 : 高斯消元·一

时间限制: 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：当然可以了，这样的小问题怎么能难到我呢？

   

提示：高斯消元

 
输入

第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

小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;
}


leibniz_zhang
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
hihocoder和leetcode-hiho-and-leet:hiho-and-leet
07-06
hihocoder 和 leetcode ##关于 Leetcode 和 Hihocoder 的一些问题和解决方案
hiho一下 第五十六周
zxiaopp的专栏
08-01 425
 题目1 : 高斯消元·一 时间限制:10000ms 单点时限:1000ms 内存限制:256MB 描述 小Ho：喂不得了啦，那边便利店的薯片半价了! 小Hi：啥?! 小Ho：那边的便利店在打折促销啊。 小Hi：走走走，赶紧去看看=v= 于是小Hi和小Ho来到了便利店。 老板为了促销，推出了组合包的形式，将不同数量的各类商品打包
模板-高斯消元(hihocoder 1195)
A Cube Master
08-05 497
高斯消元裸题，存为模板吧。值得注意的是精度1e-8、1e-10过不去，但是1e-6能过，据说每次找最大的行消元可以减小误差。 #include #include #include #include #include #include #include #include using namespace std; const double eps = 1e-6; /
hihocoder 1195&1196
Heret1c
05-18 287
1195 1196存一发高斯消元模板。（填坑） 1195是纯模板，1196是解异或方程组。 代码中有注释。【代码】 1195 #include <cstdio> #include <iostream> #include <algorithm> #include <cstring> #include <queue> #include <cmath> #define N 505 #define M
Codeforces 1195F Geometers Anonymous Club
Gipsy
07-18 480
可以写树状数组，这里直接写的线段树 #include<bits/stdc++.h> using namespace std; const int N = 1e5 + 10; struct SegmentTree { int l, r, sum, f; #define l(x) Tree[x].l #define r(x) Tree[x].r #def...
高斯消元 hihoCoder1195 高斯消元·一
qwb的博客
09-08 1203
传送门：点击打开链接 思路：高斯消元。。 没学线代前看了好久的高斯消元都没看懂，感觉线代书里的高斯消元讲的详细多了。。 然后网上找过关于高斯消元的代码，几乎没有能返回是否有解，有几组解等的代码-_- 然后今天写了一天，也总算完成了 大概的思路也就是做出梯形矩阵，然后判断矩阵的秩，详细推荐参考线性代数的书会比较好 如果不想弄懂也行，，反正模板已弄好了，以后可以直接套模板 下面是模板in
hihocoder #1195 : 高斯消元·一
yeqm_的博客
07-13 380
#include<fstream> #include<iostream> #include<string.h> #include<math.h> using namespace std; #define FOR(i,b,e) for(int i=(b);i<=(e);i++) #define FORE(i,b,e) for(int i=(b);i>=(e);i--) #define m
matlab代码做游戏-PAR-11-vs-HiHo-9-Scoring-for-Squash:蒙特卡洛模拟，观察统计数据，以将点对点计分系统
05-21
vs.HiHo-9得分 蒙特卡洛模拟，观察统计数据，以将点对点计分系统与MA Squash在2009年提出的上手交接式计分系统进行对比。 以下代码在计算赢得一场比赛的概率时假设一个简单的参数分布。 该剧本在两个玩家之间进行...
hiho-pytorch-base
02-17
hiho-pytorch-base
leetcode答案-CodingProblemReview:这是一个个人leetcode，hiho，kikstart问题解决和文书工作库，
06-30
leetcode,hiho, kikstart problem solving and paper work repository After work like a year, I am thinking my objective is cloud computing, focusing on storage, buildout, deployment. It is a little bit ...
leetcode和codewar-code-practice:肥超的刷题记录/代码锻炼
07-07
之前注意先处理一下类名和包名。 src/algorithms 文件夹内放了一些常见算法/数据结构的实现代码（目前还比较少）。 src/leetcode 文件夹内放了 leetcode 的刷题记录。 src/codewars 文件夹放了 codewars 的刷题记录...
一个圆上的n个不同点可以连成多少条不相交的线段
leibniz_zhang的博客
07-24 5129
#include #include #include using namespace std; #define LL __int64 LL pp(LL xx) { if (xx==1) return 0; if (xx==2) return 0; if (xx==3) return 0; if (xx==4)
2016年中国大学生程序设计竞赛（合肥）-重现赛（感谢安徽大学） A E H I
leibniz_zhang的博客
11-05 1723
5931 传递 Time Limit: 12000/6000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others) Total Submission(s): 71    Accepted Submission(s): 33 Problem Description 我们称一个有向图G是传递的，当且仅当对任意三个不
51nod 1165 整边直角三角形的数量 【数学：公式--求约数】
leibniz_zhang的博客
12-25 1281
1165 整边直角三角形的数量 基准时间限制：2 秒 空间限制：131072 KB 分值: 160 难度：6级算法题  收藏  关注 直角三角形，三条边的长度都是整数。给出周长N，求符合条件的三角形数量。 例如：N = 120，共有3种不同的满足条件的直角3角行。分别是：{20,48,52}, {24,45,51}, {30,40,50}。
【zzulioj】河南多校训练(985专场）<位运算--DP---思维--扩展欧几里德---回溯---字典树---状态优化>
leibniz_zhang的博客
08-03 1234
Problem B: 985的数学难题 Time Limit: 2 Sec  Memory Limit: 128 MB Submit: 85  Solved: 13 SubmitStatusWeb Board Description 985有n个正整数，他想快速知道下面函数的返回值 int a[N+1]; long long Solve() {     int i,
河南多校--可怜张继科没夺冠【打表--素数筛法--搜索--数学---物理】
leibniz_zhang的博客
08-13 1170
PS: 一些小牢骚--- 想搞好ACM--- 你不仅要会---数据结构--图论--DP---数论---各种定理---- 你还要会各种----DP优化----状态压缩------ 你更要精通---数学---物理---英语--- -----要上知天文---下知地理---- 要不然，给你个题，你都不知道是让干嘛呢--- 啊---好烦好烦-.-我要出去玩-.-
51nod oj 1137 1242 <结构体内重载运算符求矩阵乘法>
leibniz_zhang的博客
08-07 1049
链接：1137 重载运算符--重新定义运算符的作用--使返回自己想要的结果-.- 代码： #include #include #include using namespace std; struct node{ int n; int a[102][102]; void clear(){ n=0; memset(a,0,s
河南多校大一训练赛 <思维--DP--背包--最短路--最大生成树>
leibniz_zhang的博客
07-30 982
每次比赛开始的状态都好差----开始的思路还不对-- 无奈--被虐成皮了---- A - F是签到题 Time Limit:2000MS     Memory Limit:32768KB     64bit IO Format:%lld & %llu Submit Status Practice LightOJ 1198 uDebug Descript
HYSBZ 2190 hdoj 2841 2824 《（打表）欧拉和求素数对+队列求1-n中与k互质个数》
leibniz_zhang的博客
07-26 947
仪仗队 Time Limit: 10000MS   Memory Limit: 265216KB   64bit IO Format: %lld & %llu Submit Status Description 　　作为体育委员，C君负责这次运动会仪仗队的训练。仪仗队是由学生组成的N * N的方阵，为了保证队伍在行进中整齐划一，C

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值