Gauss

最新推荐文章于 2025-07-26 17:53:04 发布

转载最新推荐文章于 2025-07-26 17:53:04 发布 · 91 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/Rlemon/archive/2013/04/28/3048344.html

本文详细介绍了高斯消元法及其改进版本——高斯-若尔当消元法的实现原理与应用案例，特别关注如何高效解决含有自由变量的线性方程组问题。

  1 //Guass_未测试版本
  2 /*
  3 如果可以肯定有唯一解,可以用化成上三角，时间有一点优势;但是如果存在自由变量
  4 那么就不能这么写，而且化成上三角后，对应的行号和列号会发生错位；
  5  
  6 */
  7 const double eps=1e-10;
  8 typedef double Matrix[N][N];
  9 void gauss(Matrix A,int n){ //n个方程组，n个自由变量；A[i][n]放常量y[i]； 
 10     int i,j,r,c,id;
 11     for (r=0;r<n;r++){
 12         id=r;
 13         for (i=r+1;i<n;i++) if (fabs(A[id][r])<fabs(A[i][r])) id=i;
 14         if (id!=r) for (j=r;j<=n;j++) swap(A[id][j],A[r][j]);
 15         
 16         for (i=r+1;i<n;i++){
 17             double tmp=A[i][r]/A[r][r];
 18             for (j=r;j<=n;j++){
 19                 A[i][j]-=tmp*A[r][j];
 20             } 
 21         }
 22     }
 23     //在A[i][n]放x[i]; 
 24     for (i=n-1;i>=0;i--){
 25         for (j=n-1;j>i;j--){
 26             A[i][n]-=A[i][j]*A[j][n];
 27         }
 28         A[i][n]/=A[i][i];
 29     }    
 30 }
 31 /*
 32 化成斜线方程，可以避免行号和列号发生错位，但时间会稍慢一点，
 33 而且当方程组有解时，唯一解：x[i]=A[i][n]/A[i][i]
 34 无穷多解，A[i][n]/A[i][i]表示在自由变量取0时的x[i]; 
 35 */ 
 36 int gauss_jordan(Matrix A,int n){
 37     int i,j,r,c,id;
 38     for (r=0;r<n;r++){
 39         id=r;
 40         for (i=r+1;i<n;i++) if (fabs(A[id][r])<fabs(A[i][r])) id=i;
 41         if (fabs(A[id][r])<eps) continue;
 42         
 43         for (i=0;i<n;i++){
 44             if (i==r) continue;
 45             double f=A[i][r]/A[r][r];
 46             for (j=r+1;j<=n;j++){
 47                 A[i][j]-=f*A[i][j];
 48             }
 49         }
 50     }
 51     int cnt=0;
 52     for (i=0;i<n;i++){
 53         if (fabs(A[i][i])<eps && fabs(A[i][n])>eps) return -1;//无解；
 54         if (fabs(A[i][i])<eps && fabs(A[i][n])<eps) cnt++;//自由变量； 
 55     }    
 56     if (cnt==0) return 0;//唯一解； 
 57     return cnt;//自由元的数量； 
 58 } 
 59 //提高精度版本；行列变换用它们的最小公倍数，所以只在最后求解的时候除法；
 60 //注意数据会不会越界；
 61 int gauss_jordan(Matrix A,int n){
 62     int i,j,r,c,id;
 63     for (r=0;r<n;r++){
 64         id=r;
 65         for (i=r+1;i<n;i++) if (abs(A[id][r])<abs(A[i][r])) id=i;
 66         if (abs(A[id][r])==0) continue;
 67         
 68         for (i=0;i<n;i++){
 69             if (i==r) continue;
 70             int g=__gcd(A[r][r],A[i][r]);
 71             int tp1=A[r][r]/g,tp2=A[i][r]/g;
 72             for (j=r;j<=n;j++){
 73                 A[i][j]=A[i][j]*tp1-A[r][j]*tp2;
 74             }
 75         }
 76     }
 77     int cnt=0;
 78     for (i=0;i<n;i++){
 79         if (abs(A[i][i])==0 && abs(A[i][n])>0) return -1;//无解；
 80         if (abs(A[i][i])==0 && abs(A[i][n])==0) cnt++;//自由变量； 
 81     }    
 82     if (cnt==0) return 0;//唯一解； 
 83     return cnt;//自由元的数量； 
 84 } 
 85 
 86 
 87 /*二进制的gauss*/
 88 int gauss(Matrix A,int n,int m){
 89     int i,j,r,c,id;
 90     for (r=0,c=0;r<n && c<m;r++,c++){
 91         
 92         for (i=r;i<n;i++) if (A[i][c]) { id=i;break; }
 93         if (A[id][c])==0){ r--;continue; }
 94         if (id!=r) for (j=c;j<=m;j++) swap(A[id][j],A[r][j]);
 95         
 96         for (i=r+1;i<n;i++)if (A[i][c]){
 97             for (j=c;j<=m;j++) A[i][j]^=A[r][j];
 98         }    
 99     }
100     for (i=r;i<n;i++) if (abs(A[i][n])>0) return -1;//无解；
101     if (m-r>0) return m-r;//自由元的个数；
102      
103 }

hust 1651 **

View Code

  1 /*
  2 题意：n*m的矩阵，要么是非负整数，要么是-1，-1表示要填的数，
  3 然后告诉你每一行，每一列的和，问你有无解，每行每列至多2个-1；
  4 
  5 很简单的想法就是列出n+m个方程组然后gauss,但gauss给我们的印象就是o(n^3;
  6 所以也许就把这个想法抛弃了，但实际是每个未知数最多出现在两个方程里,也就是
  7 在行里变换中，只需要变换一行后就可以直接break，时间复杂度就是O(n^2);
  8 
  9 */
 10 #include<iostream>
 11 #include<cmath>
 12 #include<vector>
 13 #include<algorithm>
 14 #include<cstdio>
 15 #include<cstdlib>
 16 #include<cstring>
 17 using namespace std;
 18 const int N=1000+10;
 19 typedef double Matrix[N][N];
 20 int mz[N][N],c[N][N];
 21 int a[N],b[N];
 22 int n,m,sz;
 23 const double eps=1e-10;
 24 Matrix A;
 25 int gauss(Matrix A,int n,int m){
 26     int i,j,r,c,id;
 27     for (r=0,c=0;r<n && c<m;r++,c++){
 28         id=r;
 29         for (i=r;i<n;i++) if (fabs(A[i][c])>fabs(A[id][c])) id=i;
 30         if (fabs(A[id][c])<eps) {  r--;continue; }
 31         if (id!=r) for (j=c;j<=m;j++) swap(A[id][j],A[r][j]);
 32 
 33         for (i=r+1;i<n;i++){
 34             if (A[i][c]==0) continue;
 35             double f=A[i][c]/A[r][c];
 36             for (j=c;j<=m;j++){
 37                 A[i][j]-=f*A[r][j];
 38             }
 39             break;
 40         }
 41     }
 42     for (i=r;i<n;i++) if (abs(A[i][m])>0) return -1;
 43     if (m-r>0) return 2;
 44     return 1;
 45 }
 46 void work(){
 47     memset(A,0,sizeof(A));
 48     for (int i=0;i<n;i++){
 49         int sum=0;
 50         for (int j=0;j<m;j++){
 51             if (c[i][j]!=-1) A[i][c[i][j]]=1;
 52             sum+=mz[i][j];
 53         }
 54         A[i][sz]=a[i]-sum;
 55     }
 56     for (int j=0;j<m;j++){
 57         int sum=0;
 58         for (int i=0;i<n;i++){
 59             if (c[i][j]!=-1) A[j+n][c[i][j]]=1;
 60             sum+=mz[i][j];
 61         }
 62         A[n+j][sz]=b[j]-sum;
 63     }
 64   /*  cout<<"++++ "<<endl;
 65     for (int i=0;i<n+m;i++){
 66         for (int j=0;j<=sz;j++){
 67             cout<<A[i][j]<<" ";
 68         }cout<<endl;
 69     }
 70     cout<<"+++++ "<<endl;
 71     */int t=gauss(A,n+m,sz);
 72     if (t==-1) printf("No solution\n");
 73     else if (t==1) printf("Unique\n");
 74     else printf("More than one\n");
 75 }
 76 int main(){
 77     int T,cas=0;
 78     scanf("%d",&T);
 79     while (T--){
 80         scanf("%d%d",&n,&m);
 81         sz=0;
 82         memset(c,-1,sizeof(c));
 83         for (int i=0;i<n;i++){
 84             for (int j=0;j<m;j++){
 85                 scanf("%d",&mz[i][j]);
 86                 if (mz[i][j]==-1){
 87                     mz[i][j]=0;
 88                     c[i][j]=sz++;
 89                 }
 90             }
 91         }
 92         for (int i=0;i<n;i++) scanf("%d",&a[i]);
 93         for (int j=0;j<m;j++) scanf("%d",&b[j]);
 94         printf("Case #%d: ",++cas);
 95         work();
 96 
 97 
 98     }
 99 
100     return 0;
101 }
102 /*
103 5
104 2 2
105 1 1
106 -1 -1
107 2 4
108 3 3
109 2 2
110 1 1
111 -1 -1
112 2 5
113 3 3
114 2 2
115 1 1
116 2 2
117 2 4
118 3 3
119 4 4
120 1 2 3 4
121 1 -1 -1 2
122 3 -1 -1 4
123 5 6 7 8
124 10 6 14 26
125 10 13 15 18
126 2 2
127 1 -1
128 2 -1
129 -1 1
130 3 -3
131 
132 */

转载于:https://www.cnblogs.com/Rlemon/archive/2013/04/28/3048344.html