首页 > 代码库 > 高斯消元 模板

高斯消元 模板

照着czyuan的那个模板,手敲了一遍,存一下。

貌似今天一整天就看了一下高斯消元的知识,然后看了模板,又手敲了一遍。

  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         else126         {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 }