您的位置:首页 > 编程语言 > C语言/C++

高斯消元模板

2015-05-21 08:39 204 查看
#include <iostream>

#include <cstdio>

#include <cstring>

/*

高斯消元就是模拟矩阵的线性变换,把矩阵化成上三角 然后判断解的情况.

*/

typedef __int64 ll;

using namespace std;

const int maxn = 105;

ll a[maxn][maxn],x[maxn];//增广矩阵

int row,col;//线性方程组的方程数为row 变元数为col 则增广矩阵的行数为 row(0-(row-1)) col(0-col)

bool free_x[maxn];

ll gcd(ll x,ll y)

{

    if(!y) return x;

    return gcd(y,x%y);

}

ll lcm(ll x,ll y)

{

    return x/gcd(x,y)*y;

}

void Debug()

{

    for(int i=0;i<row;i++)

    {

        for(int j=0;j<col+1;j++)

            printf("%I64d ",a[i][j]);

        printf("\n");

    }

    printf("\n");

}

void swaprowxy(int x,int y)//交换x y两行

{

    for(int i=0;i<col+1;i++)

    {

        int tp=a[x][i];

        a[x][i]=a[y][i];

        a[y][i]=tp;

    }

}

void swapcolxy(int x,int y)//交换x y两列

{

    for(int i=0;i<row;i++)

    {

        int tp=a[i][x];

        a[i][x]=a[i][y];

        a[i][y]=tp;

    }

}

ll Abs(ll x)

{

    if(x<0) return -x;

    return x;

}

int Gauss()

{

    int r,c;

    for(r=0,c=0;r<row&&c<col;r++,c++) //用第r行消掉其余((r+1)-(row-1))行的第c列

    {

        int maxr=r;//首先找到 (r-(row-1))行中第c列绝对值最大的行 与第r行交换(防止除法误差什么的?)

        for(int i=r+1;i<row;i++)

            if(Abs(a[i][c])>Abs(a[maxr][c])) maxr=i;

        if(maxr!=r)

            swaprowxy(maxr,r);

        //当把第r行换成Abs(a[r][c])最大的时候仍然出现 a[r][c]==0 则说明以下的所有行的第c列都为0 则直接处理第r行的第c+1列

        if(a[r][c]==0)

        {

            r--;

            continue;

        }

        for(int i=r+1;i<row;i++)//枚举以下的所有行 利用第r行消除掉以下所有行的第c列

        {

            ll lc=lcm(Abs(a[r][c]),Abs(a[i][c]));

            ll t1=lc/a[r][c],t2=lc/a[i][c];

            int sign=1;

            if((a[r][c]<0&&a[i][c]>0)||(a[r][c]>0&&a[i][c]<0))

                sign=-1;

            for(int j=c;j<col+1;j++)

            {

                if(!a[i][j]) continue;

                a[i][j]=a[i][j]*t2-sign*a[r][j]*t1;

            }

        }

    }

    Debug();

    //1.return -1; 这种情况下是无解的 说明增广矩阵中存在着(0,0....,a)这样的行(a!=0)

    for(int i=r;i<row;i++)

        if(a[i][c]!=0) return -1;

    //2.return -2; 这种情况下有无穷解 在增广矩阵中出现(0,0...,0)这样的行 没有形成严格的上三角

    //出现的行数就等于自由变元的个数

    if(r<col)

    {

        for(int i=r-1;i>=0;i--)

        {

            int free_num=0,free_index;//变元的数量和变元的下标

            for(int j=0;j<col;j++)

            {

                if(a[i][j]!=0&&!free_x[j]) free_num++,free_index=j;

            }

            if(free_num>1) continue;//未确定的变元数量大于1

            //这个方程中只有一个x是未知的

            int temp=a[i][col];//第i个方程的常数

            for(int j=0;j<col;j++)

            {

                if(a[i][j]!=0&&j!=free_index) temp-=a[i][j]*x[j];

            }

            x[free_index]=temp/a[i][free_index];

            free_x[free_index]=1;//这个x的值是确定的

        }

        return col-r;

    }

    for(int i=col-1;i>=0;i--)

    {

        int temp=a[i][col];

        for(int j=i+1;j<col;j++)

        {

            if(a[i][j]!=0) temp-=a[i][j]*x[j];

        }

        if(temp%a[i][i]!=0) return -2;//浮点数解 没有整数解.

        x[i]=temp/a[i][i];

    }

    return 0;

}

int main()

{

    //freopen("in.txt","r",stdin);

    while(~scanf("%d%d",&row,&col))

    {

        memset(a,0,sizeof(a));

        memset(free_x,0,sizeof(free_x)); //初始所有的元都是不确定的

        memset(x,0,sizeof(x));//记录xi的值

        for(int i=0;i<row;i++)

            for(int j=0;j<col+1;j++)

                scanf("%I64d",&a[i][j]);

        Debug();

        int ans_num=Gauss(); //返回解的个数

        printf("%d\n",ans_num);

    }

    return 0;

}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签:  C++ 算法