高斯消元 模板

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <cstdlib>
  5 #include <cmath>
  6 #include <algorithm>
  7 #define LL __int64
  8 const int maxn = 100+10;
  9 using namespace std;
 10 int equ, var, fn;
 11 int a[maxn][maxn], x[maxn];
 12 bool free_x[maxn];
 13 
 14 int gcd(int a, int b)
 15 {
 16     return b==0?a:gcd(b, a%b);
 17 }
 18 int lcm(int a, int b)
 19 {
 20     return a*b/gcd(a, b);
 21 }
 22 int Gauss()
 23 {
 24     int i, j, k, max_r, col;
 25     int ta, tb, LCM, tmp, fx_num;
 26     int free_index;
 27     col = 0;
 28 
 29     for(k = 0; k<equ && col<var; k++, col++)
 30     {
 31         max_r = k;
 32         for(i = k+1; i < equ; i++)
 33             if(abs(a[i][col])>abs(a[max_r][col]))
 34             max_r = i;
 35 
 36         if(max_r != k)
 37             for(j = k; j < var+1; j++)
 38             swap(a[k][j], a[max_r][j]);
 39 
 40         if(a[k][col]==0)
 41         {
 42             k--;
 43             continue;
 44         }
 45         for(i = k+1; i < equ; i++)
 46         {
 47             if(a[i][col] != 0)
 48             {
 49                 LCM = lcm(abs(a[i][col]), abs(a[k][col]));
 50                 ta = LCM/abs(a[i][col]);
 51                 tb= LCM/abs(a[k][col]);
 52                 if(a[i][col]*a[k][col] < 0) tb = -tb;
 53 
 54                 for(j = col; j < var+1; j++)
 55                 a[i][j] = a[i][j]*ta - a[k][j]*tb;
 56             }
 57         }
 58     }
 59     for(i = k; i < equ; i++)
 60     if(a[i][col] != 0)
 61     return -1;
 62 
 63     if(k < var)
 64     {
 65         for(i = k-1; i >= 0; i--)
 66         {
 67             fx_num = 0;
 68             for(j = 0; j < var; j++)
 69                 if(a[i][j]!=0 && free_x[j])
 70                 {
 71                     fx_num ++;
 72                     free_index = j;
 73                 }
 74             if(fx_num > 1) continue;
 75 
 76             tmp = a[i][var];
 77             for(j = 0; j < var; j++)
 78                 if(a[i][j] != 0 && j!=free_index)
 79                 tmp -= a[i][j]*x[j];
 80 
 81             x[free_index] = tmp/a[i][free_index];
 82             free_x[free_index] = 0;
 83         }
 84         return var-k;
 85     }
 86     for(i = var-1; i >= 0; i--)
 87     {
 88         tmp = a[i][var];
 89         for(j = i+1; j < var; j++)
 90         if(a[i][j] != 0)
 91         tmp -= a[i][j]*x[j];
 92 
 93         if(tmp%a[i][i] != 0) return -2;
 94         x[i] = tmp/a[i][i];
 95     }
 96     return 0;
 97 }
 98 int main()
 99 {
100     int i, j;
101     while(cin>>equ>>var)
102     {
103         memset(a, 0, sizeof(a));
104         memset(x, 0, sizeof(x));
105         memset(free_x, 1, sizeof(free_x));
106         for(i = 0; i < equ; i++)
107         {
108             for(j = 0; j < var+1; j++)
109             cin>>a[i][j];
110         }
111 
112         fn = Gauss();
113         if(fn == -1) cout<<"无解"<<endl;
114         else if(fn == -2)
115         cout<<"有浮点数解， 无整数解"<<endl;
116         else if(fn > 0)
117         {
118             printf("无穷多解! 自由变元个数为%d\n", fn);
119             for(i = 0; i < var; i++)
120             {
121                 if(free_x[i]) printf("x%d行 是不确定的\n", i + 1);
122                 else printf("x%d: %d\n", i + 1, x[i]);
123             }
124         }
125         else
126         {
127             for(i = 0; i < var; i++)
128             printf("x%d: %d\n", i + 1, x[i]);
129         }
130         cout<<endl;
131     }
132     return 0;
133 }