高斯消元法模板

来源:互联网 发布:2017网络新技术 编辑:程序博客网 时间:2024/06/04 19:35

浮点数版(以hdu 3359为例):

#pragma comment(linker, "/STACK:102400000,102400000")#include <stdio.h>#include <iostream>#include <algorithm>#include <sstream>#include <stdlib.h>#include <string.h>#include <limits.h>#include <string>#include <time.h>#include <math.h>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define INF 0x3f3f3f3f#define eps 1e-8#define pi acos(-1.0)typedef long long ll;double s[20][20],a[200][200],x[200];int equ,var;int Guass(){//初始化时等号右边的常数放到x数组里,计算结果保存在x数组里。int i,j,k,col,max_r;for(k=0,col=0;k<equ&&col<var;k++,col++){max_r=k;for(i=k+1;i<equ;i++)if(fabs(a[i][col])>fabs(a[max_r][col]))max_r=i;if(fabs(a[max_r][col])<eps)return 0;if(k!=max_r){for(j=col;j<var;j++)swap(a[k][j],a[max_r][j]);swap(x[k],x[max_r]);}x[k]/=a[k][col];for(j=col+1;j<var;j++)a[k][j]/=a[k][col];a[k][col]=1;for(i=0;i<equ;i++)if(i!=k){x[i]-=x[k]*a[i][k];for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col];a[i][col]=0;}}return 1;}int main(){     //freopen("data.in","r",stdin);     //freopen("data.out","w",stdout); int d,tt=0,m,n;     while(~scanf("%d%d%d",&m,&n,&d)&&(n+m+d)){ if(tt++)puts(""); for(int i=0;i<n;i++) for(int j=0;j<m;j++) scanf("%lf",&s[i][j]); memset(a,0,sizeof(a)); equ=var=n*m; for(int i=0;i<n;i++) for(int j=0;j<m;j++){ int cnt=0; for(int k=0;k<n;k++) for(int p=0;p<m;p++) if(abs(i-k)+abs(j-p)<=d){ a[i*m+j][k*m+p]=1; cnt++; } x[i*m+j]=cnt*s[i][j]; } Guass();// cout<<ha<<endl; for(int i=0;i<n;i++){ for(int j=0;j<m;j++)  printf ("%8.2lf", x[i*m+j]); puts(""); } }     return 0;}


整数版:

/* ***********************************************Author :rabbitCreated Time :2014/3/7 15:19:27File Name :11.cpp************************************************ */#pragma comment(linker, "/STACK:102400000,102400000")#include <stdio.h>#include <iostream>#include <algorithm>#include <sstream>#include <stdlib.h>#include <string.h>#include <limits.h>#include <string>#include <time.h>#include <math.h>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define INF 0x3f3f3f3f#define eps 1e-8#define pi acos(-1.0)typedef long long ll;const int maxn=50;int a[maxn][maxn];//增广矩阵。int x[maxn];//解集。bool free_x[maxn];//标记是否为不确定的变元。int gcd(int a,int b){if(a==0)return b;return gcd(b%a,a);}int lcm(int a,int b){return a/gcd(a,b)*b;}int Guass(int equ,int var){int max_r;//当前这一列最大的行。int ta,tb,LCM,temp,free_x_num,free_index,i,j,k,col;for(i=0;i<=var;i++)x[i]=0,free_x[i]=true;for(col=0,k=0;k<equ&&col<var;k++,col++){max_r=k;for(i=k+1;i<equ;i++)if(abs(a[i][col])>abs(a[max_r][col]))max_r=i;if(max_r!=k){for(j=k;j<var+1;j++)swap(a[k][j],a[max_r][j]);}if(a[k][col]==0){k--;continue;//col列第k行之下全是0.处理当前行的下一列。}for(i=k+1;i<equ;i++)if(a[i][col]){LCM=lcm(abs(a[i][col]),abs(a[k][col]));ta=LCM/abs(a[i][col]);tb=LCM/abs(a[k][col]);if(a[i][col]*a[k][col]<0)tb=-tb;for(j=col;j<var+1;j++)a[i][j]=a[i][j]*ta-a[k][j]*tb;}}for(i=k;i<equ;i++)if(a[i][col])return -1;//无解。if(k<var){for(i=k-1;i>=0;i--){free_x_num=0;for(j=0;j<var;j++)if(a[i][j]&&free_x[j])free_x_num++,free_index=j;if(free_x_num>1)continue;//无法求解出确定变元。temp=a[i][var];for(j=0;j<var;j++)if(a[i][j]&&j!=free_index)temp-=a[i][j]*x[j];x[free_index]=temp/a[i][free_index];//求出变元。free_x[free_index]=0;//该变元是确定的。}return var-k;//自由变元有var-k个。}for(i=var-1;i>=0;i--){temp=a[i][var];for(j=i+1;j<var;j++)if(a[i][j])temp-=a[i][j]*x[j];if(temp%a[i][i])return -2;//有浮点数解,但是无整数解。x[i]=temp/a[i][i];}return 0;}int main(){    int i, j;    int equ,var;    while (scanf("%d %d", &equ, &var) != EOF){        memset(a, 0, sizeof(a));        for (i = 0; i < equ; i++)            for (j = 0; j < var + 1; j++)                scanf("%d", &a[i][j]);        int free_num = Guass(equ,var);        if (free_num == -1) printf("无解!\n");        else if (free_num == -2) printf("有浮点数解,无整数解!\n");        else if (free_num > 0){            printf("无穷多解! 自由变元个数为%d\n", free_num);            for (i = 0; i < var; i++){                if (free_x[i]) printf("x%d 是不确定的\n", i + 1);                else printf("x%d: %d\n", i + 1, x[i]);            }        }        else{            for (i = 0; i < var; i++)                printf("x%d: %d\n", i + 1, x[i]);        }        printf("\n");    }    return 0;}/*4 42 1 -5 1 81 -3 0 -6 90 2 -1 2 -51 4 -7 6 04 41 1 1 1 51 2 -1 4 -22 -3 -1 -5 -23 1 2 11 0*/ 

位运算高斯消元(以poj  1222为例):

代码:

/* ***********************************************Author :xianxingwuguanCreated Time :2014/3/7 14:25:12File Name :11.cpp************************************************ */#pragma comment(linker, "/STACK:102400000,102400000")#include <stdio.h>#include <iostream>#include <algorithm>#include <sstream>#include <stdlib.h>#include <string.h>#include <limits.h>#include <string>#include <time.h>#include <math.h>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define INF 0x3f3f3f3f#define eps 1e-8#define pi acos(-1.0)typedef long long ll;const int maxn=50;int a[maxn][maxn],x[maxn],free_num,free_x[maxn];int Guass(int equ,int var){//-1返回无解,0返回唯一解,否则返回自由变元个数,int i,j,col,k,free_num=0;memset(x,0,sizeof(x));for(k=0,col=0;k<equ&&col<var;k++,col++){int max_r=k;for(i=k+1;i<equ;i++)if(abs(a[i][col])>abs(a[max_r][col]))max_r=i;if(max_r!=k){for(j=k;j<var+1;j++)swap(a[k][j],a[max_r][j]);}if(a[k][col]==0){k--;free_x[free_num++]=col;//自由变元continue;}for(i=k+1;i<equ;i++)if(a[i][col]){for(j=col;j<var+1;j++)a[i][j]^=a[k][j];}}for(i=k;i<equ;i++)if(a[i][col])return -1;//无解if(k<var)return var-k;//自由变元个数for(int i=var-1;i>=0;i--){x[i]=a[i][var];for(int j=i+1;j<var;j++)x[i]^=(a[i][j]&x[j]);}return 0;}int main(){     //freopen("data.in","r",stdin);     //freopen("data.out","w",stdout);  int T,u,v,n; scanf("%d",&T); for(int t=1;t<=T;t++){ printf("PUZZLE #%d\n",t); memset(a,0,sizeof(a)); int var=30,equ=30; for(int i=0;i<5;i++) for(int j=0;j<6;j++){ int p=i*6+j; a[p][p]=1; if(i>0)a[(i-1)*6+j][p]=1; if(i<4)a[(i+1)*6+j][p]=1; if(j>0)a[i*6+j-1][p]=1; if(j<5)a[i*6+j+1][p]=1; }         for(int i = 0;i < 30;i++)scanf("%d",&a[i][30]);         Guass(equ,var);         for(int i = 0;i < 5;i++) {             for(int j = 0;j < 5;j++)printf("%d ",x[i*6+j]);             printf("%d\n",x[i*6+5]);         } }     return 0;}


0 0