高斯消元

来源：互联网发布：周扬青淘宝编辑：程序博客网时间：2024/06/16 20:57
手打高斯未优化以后补充
解释以后补充
#include<iostream>#include<cstdio>#include<iomanip>#include <cmath>using namespace std;int d[50][50];int n,m;int gcd(int a,int b){    return a%b==0?b:gcd(b,a%b);}void rungao(){    for(int i=0;i<n;i++)    {        for(int j=i+1;j<m;j++)        {            if(!d[i][j]) continue;            int k=j;                int x=(d[i][j]);                int y=(d[k][k]);                if(x&&y)                {                    int z=gcd(x,y);                    int x1=x*y/z/x,y1=x*y/z/y;                    //cout<<x1<<' '<<y1<<endl;                    int re=0;                    for(int u=0;u<m;u++)                    {                        d[i][u]=x1*d[i][u]-y1*d[k][u];                        if(d[i][u])                        {                            if(!re) re=abs(d[i][u]);                            re=gcd(re,abs(d[i][u]));                        }                    }                    if(!re) re=!re;                    for(int u=0;u<m;u++)                    {                        d[i][u]/=re;                        //cout<<d[i][u]<<' ';                    }//cout<<endl;                }        }    }}void gaosi(){    for(int i=0,j=0;i<n&&j<m;i++,j++)    {        if(!d[i][j]) swap(d[i],d[i+1]);        for(int k=i+1;k<n;k++)        {            int x=d[i][j];            int y=d[k][j];            if(x&&y)            {                int z=gcd(x,y);                int x1=x*y/z/y,y1=x*y/z/x;                //cout<<x1<<' '<<y1<<endl;                for(int u=j;u<m;u++)                {                    d[k][u]=x1*d[k][u]-y1*d[i][u];                    //cout<<d[k][u]<<' ';                }                //cout<<endl;                int re=0;                for(int u=j;u<m;u++)                {                   // d[k][u]=x1*d[k][u]-y1*d[i][u];                    if(d[k][u])                    {                        if(re==0) re=abs(d[k][u]);                        //cout<<re<<' ';                        re=gcd(re,abs(d[k][u]));                    }                    //cout<<d[k][u]<<' ';                }                if(!re) re=!re;                //cout<<endl;                for(int u=j;u<m;u++)                {                    d[k][u]/=re;                    //cout<<d[k][u]<<' ';                }                //cout<<endl;            }        }    }        for(int i=0;i<n;i++)        {            for(int j=0;j<m;j++)            {                cout<<setw(7)<<d[i][j]<<' ';            }            cout<<endl;        }        cout<<endl;    rungao();}int main(){    while(cin>>n>>m)    {        for(int i=0;i<n;i++)        {            for(int j=0;j<m;j++)            {                cin>>d[i][j];            }        }gaosi();        cout<<endl;        for(int i=0;i<n;i++)        {            for(int j=0;j<m;j++)            {                cout<<setw(7)<<d[i][j]<<' ';            }            cout<<endl;        }    }}
阅读全文
0 0