现在的位置: 首页 > 综合 > 正文

高斯消元

2013年03月27日 ⁄ 综合 ⁄ 共 5676字 ⁄ 字号 小 中 大 ⁄ 评论关闭

高斯消元法（理论+习题）

高斯消元主要用来解多元一次线性方程组，主要用到的线性代数的知识详见矩阵行列式那一篇       （http://blog.csdn.net/zhoufenqin/article/details/7779707）

求解线性方程的结果会出现三种情况：无解，多解和唯一解。看下图

        

当某一行出现（0,0，……0，a) 时，方程无解。因为x1*0+x2*0+……+xn*0=a（a不为0）等式无解

当从某一行开始至最后一行，出现（0,0，……0）说明方程存在自由元，则会出现多解

当矩阵的秩等于等式个数时，有唯一解


        举两个关于高斯消元的问题：
       一.经典的画家问题（poj 1681Painter's Problem）
题意：在由n*n的方块组成的画板上，每个格子不是白色就是黄色，染坐标为（i，j）的格子时 (i-1, j), (i+1, j), (i, j-1) 和 (i, j+1)的格子的颜色会相应变化（原来黄色的），要把画板的所有格子都染成黄色，最少要几步
思路：
  如图所示，当改变红色方格时，黄色方格的颜色都会相应变化，同理，只要这四个黄色方格有一个发生变化时，红色方格都会变化，那么当四个黄色方格都发生变化时，红色方格最终又变成原来的颜色，所以任意一个格子最多只需要变化一次，因为变化偶数次相当于没有变化，变化奇数次的效果又跟变化一次的效果一样。
所以每个格子都只有两种情况0和1,0代表没有变化，1代表变化一次。

红色格子最后的颜色由周围四个黄色格子变化次数决定，我们假设四个黄色格子的变化次数分别为X1，X2，X3，X4，则红色格子变化=（X1+X2+X3+X4)%2;

1 2 3
4 5 6
7 8 9
从上图中可知，1格子会影响2,4两个格子，2格子会影响1,5,3；5格子会影响2,4,6,8，以此类推……因此我们可以列出以下矩阵：
第一行中为1的数代表第一个方格对1,2,4三个方格有影响，同理第二个方格对1,2,3,5方格有影响………………
于是我们列出9*9的举证来表述方格之间的相互影响。现在把得到的矩阵变成增广矩阵，那么最后一列应该填什么呢？最后一列表示画板从原来的状态变成最后的状态，例如一个格子，按题目意思，如果原先是黄色，则最后还是黄色，不用发生变化，则第一行最后一列填0，因为(X1+X2+X3)%2=0，同理，如果原先是白色，则此方格要发生变化，则最后一列填1.

进过分析，我们就可以顺利地写出代码啦！！

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>

using namespace std;

char str[20];
int paint[400][400];
int N,st;
int Find[10]={1,0,-1,0,0,-1,0,1};

void Guass()//高斯消元
{
    int i,j,row,col;
    for(row=0,col=0;row<st && col<st;row++,col++)//从第一行开始，不停地把原矩阵化成上三角
    {
        for(i=row;i<st;i++)//找到当前行以下的不为0的一行，有些也写成找当前行以下当前列值最大的一行，无影响
            if(paint[i][col]) break;//找到就跳出
        if(i==st)//如果下面每一行的当前列已经都是0，那么从当前行的后一个数继续找（回顾线性代数的知识）
        {
            row--;
            continue;
        }
        if(i!=row)//如果找到的那一行不是当前行，则进行整行交换
        {
            for(j=0;j<st+1;j++) swap(paint[row][j],paint[i][j]);
        }
        for(i=row+1;i<st;i++)
            if(paint[i][col])
            {
                for(j=col;j<st+1;j++)//把当前行一下的当前列的值化为0
                    paint[i][j]^=paint[row][j];
            }
    }
    for(i=row;i<st;i++)//查找，如果出现（0,0,……0,a）方程无解
        if(paint[i][st]!=0)
        {
            printf("inf\n");
            return ;
        }
    int maxn=0x7fffffff;
    int Num=1<<(st-row);//自由元有st-row个，则进行枚举
    int pos;
    int cnt;
    for(i=0;i<Num;i++)//枚举的过程
    {
        cnt=0;
        for(j=st-1,pos=i;j>=row;j--,pos>>=1)
        {
            paint[j][j]=pos&1;
            if(paint[j][j]) cnt++;
        }
        for(j=row-1;j>=0;j--)
        {
            int temp=0;
            for(int k=st-1;k>j;k--)
            {
                if(paint[j][k]==0) continue;
                temp^=paint[k][k];
            }
            paint[j][j]=paint[j][st]^temp;
            if(paint[j][j]) cnt++;
        }
        maxn=min(maxn,cnt);
    }
    printf("%d\n",maxn);
}
int L(int x,int y)
{
    return x*N+y;
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        memset(paint,0,sizeof(paint));
        scanf("%d",&N);
        st=N*N;
        for(int i=0;i<N;i++)
        {
            scanf("%s",str);
            for(int j=0;j<N;j++)//建立格子之间的相互联系
            {
                paint[L(i,j)][st]=(str[j]=='y'?0:1);
                paint[L(i,j)][L(i,j)]=1;
                for(int k=0;k<4;k++)
                {
                    int x=i+Find[k];
                    int y=j+Find[k+4];
                    if(x<0 || x>=N || y<0 || y>=N)
                        continue;
                    paint[L(i,j)][L(x,y)]=1;
                }
            }
        }

        Guass();
    }
    return 0;
}
二：解同余方程式的高斯消元（poj 2947  Widget Factory）
题意：N种物品，M条记录，接写来M行，每行有K，str1，str2，表述从星期str1到星期str2，做了K件物品，接下来的K个数为物品的编号。此题注意最后结果要调整到3-9之间
    思路，很容易想到高斯消元。但是是带同余方程式的高斯消元，开始建立关系的时候就要MOD 7
解此类方程式时最后求解的过程要用到扩展gcd的思想，举个例子，如果最后得到的矩阵为：
    1  1   4
    0  6   4
   则6 * X2 % 7= 4 % 7  则相当于6 * X2 + 7 * Y = 4，利用扩展gcd则可求出X2为3，则第一个方程为
X1 * 1 % 7 + 1*3 % 7 = 4 % 7 则相当于 X1 + 7 * Y = 1  得到X1=1；

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>

using namespace std;

#define MOD 7
int a[310][310];
int N,M;//M件物品，N个记录！！
char str1[5],str2[5];
int ans[310];
int solved(char s[])//额~此处可以优化
{
    if(s[0]=='M') return 1;
    if(s[0]=='T' && s[1]=='U') return 2;
    if(s[0]=='W') return 3;
    if(s[0]=='T' && s[1]=='H') return 4;
    if(s[0]=='F') return 5;
    if(s[0]=='S' && s[1]=='A') return 6;
    if(s[0]=='S' && s[1]=='U') return 7;
}
int extend_gcd(int A,int B,int &x,int &y)//扩展gcd
{
    if(B==0)
    {
        x=1,y=0;
        return A;
    }
    else
    {
        int r=extend_gcd(B,A%B,x,y);
        int t=x;
        x=y;
        y=t-A/B*y;
        return r;
    }
}
int lcm(int A,int B)//求最小公倍数
{
    int x=0,y=0;
    return A*B/extend_gcd(A,B,x,y);
}
void Guass()//高斯消元
{
    int i,j,row,col;
    for(row=0,col=0;row<N && col<M;row++,col++)
    {
        for(i=row;i<N;i++)
            if(a[i][col]) break;
        if(i==N)
        {
            row--;
            continue;
        }
        if(i!=row)
            for(j=0;j<=M;j++) swap(a[row][j],a[i][j]);
        for(i=row+1;i<N;i++)
        {
            if(a[i][col])
            {
                int LCM=lcm(a[row][col],a[i][col]);//利用最小公倍数去化上三角
                int ch1=LCM/a[row][col],ch2=LCM/a[i][col];
                for(j=col;j<=M;j++)
                    a[i][j]=((a[i][j]*ch2-a[row][j]*ch1)%MOD+MOD)%MOD;
            }

        }

    }

    for(i=row;i<N;i++)//无解的情况
    {
        if(a[i][M]!=0)
        {
            printf("Inconsistent data.\n");
            return ;
        }
    }
    if(row<M)//有自由元的情况，则会出现多种解
    {
        printf("Multiple solutions.\n");
        return ;
    }
    //当有唯一解时
    for(i=M-1;i>=0;i--)
    {
        int ch=0;
        for(j=i+1;j<M;j++)
        {
            ch=(ch+ans[j]*a[i][j]%MOD)%MOD;
        }
        int last=((a[i][M]-ch)%MOD+MOD)%MOD;
        int x=0,y=0;
        int d=extend_gcd(a[i][i],MOD,x,y);//利用扩展gcd求解
        x%=MOD;
        if(x<0) x+=MOD;
        ans[i]=last*x/d%MOD;
        if(ans[i]<3) ans[i]+=7;
    }
    for(int i=0;i<M;i++)
    {
        if(i==0)
            printf("%d",ans[i]);
        else
            printf(" %d",ans[i]);

    }
    printf("\n");
}
int main()
{

    while(scanf("%d%d",&M,&N),N||M)
    {
        memset(a,0,sizeof(a));
        memset(ans,0,sizeof(ans));
        int K;

        for(int i=0;i<N;i++)
        {
            scanf("%d%s%s",&K,str1,str2);
            int last=((solved(str2)-solved(str1)+1)%MOD+MOD)%MOD;

            for(int j=0;j<K;j++)
            {
                int v;
                scanf("%d",&v);
                a[i][v-1]=(a[i][v-1]+1)%MOD;
            }
            a[i][M]=last;
        }

        Guass();
    }
    return 0;
}


注意

1.自由元不一定都在最后，所以有时候需要进行列变化，把自由元都调到最后去，这样才方便枚举
    for(j=0;j<N;j++)
    {
        if(A[j][j]) continue;
        for(k=j;k<N;k++)
            if(A[j][k]) break;
        if(k==N) break;
        for(i=0;i<N;i++)
            swap(A[i][j],A[i][k]);
    }
2.当行和列数目太多时，要用dfs去搜索，因为(1<<x)太大无法存下来
void dfs(int u)//一开始传进来GUASS()得到的row
{
    if(u==N)
    {
        for(int i=0;i<N;i++)
            temp[i]=ans[i];

        for(int i=row-1;i>=0;i--)
        {
            int ch=A[i][N];
            for(int j=i+1;j<N;j++)
                ch=ch^(temp[j] && A[i][j]);
            temp[i]=ch;
        }
        int sum=0;
        for(int i=0;i<N;i++)
            sum+=temp[i];
        maxn=min(maxn,sum);
        return ;
    }
    ans[u]=0;
    dfs(u+1);
    ans[u]=1;
    dfs(u+1);
}
 poj上的几道高斯消元：
返回
【上篇】
【下篇】

作者:

抱歉!评论已关闭.

返回首页

Copyright © 2013-2018 学步园  保留所有权利.
软文销售 QQ客服：2265327166 点击这里给我发消息（其他合作也可洽谈）
必威体育
必威电竞