高斯消元问题

可能是目前做的这类题都比较简单,所以几乎都是套上模板就差不多了.....

这里copy一份czyuan大神的模板

  1 /* 用于求整数解得方程组. */
  2 #include <iostream>
  3 #include <string>
  4 #include <cmath>
  5 using namespace std;
  6 const int maxn = 105;
  7 int equ, var; // 有equ个方程,var个变元。增广阵行数为equ, 分别为0到equ - 1,列数为var + 1,分别为0到var.
  8 int a[maxn][maxn];
  9 int x[maxn]; // 解集.
 10 bool free_x[maxn]; // 判断是否是不确定的变元.
 11 int free_num;
 12 void Debug(void)
 13 {
 14     int i, j;
 15     for (i = 0; i < equ; i++)
 16     {
 17         for (j = 0; j < var + 1; j++)
 18         {
 19             cout << a[i][j] << " ";
 20         }
 21         cout << endl;
 22     }
 23     cout << endl;
 24 }
 25 inline int gcd(int a, int b)
 26 {
 27     int t;
 28     while (b != 0)
 29     {
 30         t = b;
 31         b = a % b;
 32         a = t;
 33     }
 34     return a;
 35 }
 36 inline int lcm(int a, int b)
 37 {
 38     return a * b / gcd(a, b);
 39 }
 40 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)
 41 int Gauss(void)
 42 {
 43     int i, j, k;
 44     int max_r; // 当前这列绝对值最大的行.
 45     int col; // 当前处理的列.
 46     int ta, tb;
 47     int LCM;
 48     int temp;
 49     int free_x_num;
 50     int free_index;
 51     // 转换为阶梯阵.
 52     col = 0; // 当前处理的列.
 53     for (k = 0; k < equ && col < var; k++, col++)
 54     { // 枚举当前处理的行.
 55         // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
 56         max_r = k;
 57         for (i = k + 1; i < equ; i++)
 58         {
 59             if (abs(a[i][col]) > abs(a[max_r][col])) max_r = i;
 60         }
 61         if (max_r != k)
 62         { // 与第k行交换.
 63             for (j = k; j < var + 1; j++) swap(a[k][j], a[max_r][j]);
 64         }
 65         if (a[k][col] == 0)
 66         { // 说明该col列第k行以下全是0了,则处理当前行的下一列.
 67             k--; continue;
 68         }
 69         for (i = k + 1; i < equ; i++)
 70         { // 枚举要删去的行.
 71             if (a[i][col] != 0)
 72             {
 73                 LCM = lcm(abs(a[i][col]), abs(a[k][col]));
 74                 ta = LCM / abs(a[i][col]), tb = LCM / abs(a[k][col]);
 75                 if (a[i][col] * a[k][col] < 0) tb = -tb; // 异号的情况是两个数相加.
 76                 for (j = col; j < var + 1; j++)
 77                 {
 78                     a[i][j] = a[i][j] * ta - a[k][j] * tb;
 79                 }
 80             }
 81         }
 82     }
 83     Debug();
 84     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
 85     for (i = k; i < equ; i++)
 86     { // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换.
 87         if (a[i][col] != 0) return -1;
 88     }
 89     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.
 90     // 且出现的行数即为自由变元的个数.
 91     if (k < var)
 92     {
 93         // 首先,自由变元有var - k个,即不确定的变元至少有var - k个.
 94         for (i = k - 1; i >= 0; i--)
 95         {
 96             // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行.
 97             // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的.
 98             free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元.
 99             for (j = 0; j < var; j++)
100             {
101                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
102             }
103             if (free_x_num > 1) continue; // 无法求解出确定的变元.
104             // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的.
105             temp = a[i][var];
106             for (j = 0; j < var; j++)
107             {
108                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
109             }
110             x[free_index] = temp / a[i][free_index]; // 求出该变元.
111             free_x[free_index] = 0; // 该变元是确定的.
112         }
113         return var - k; // 自由变元有var - k个.
114     }
115     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
116     // 计算出Xn-1, Xn-2 ... X0.
117     for (i = var - 1; i >= 0; i--)
118     {
119         temp = a[i][var];
120         for (j = i + 1; j < var; j++)
121         {
122             if (a[i][j] != 0) temp -= a[i][j] * x[j];
123         }
124         if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
125         x[i] = temp / a[i][i];
126     }
127     return 0;
128 }
129 int main(void)
130 {
131     //freopen("Input.txt", "r", stdin);
132     int i, j;
133     while (scanf("%d %d", &equ, &var) != EOF)
134     {
135         memset(a, 0, sizeof(a));
136         memset(x, 0, sizeof(x));
137         memset(free_x, 1, sizeof(free_x)); // 一开始全是不确定的变元.
138         for (i = 0; i < equ; i++)
139             for (j = 0; j < var + 1; j++)
140                 scanf("%d", &a[i][j]);
141         //        Debug();
142         free_num = Gauss();
143         if (free_num == -1) printf("无解!
");
144         else if (free_num == -2) printf("有浮点数解,无整数解!
");
145         else if (free_num > 0)
146         {
147             printf("无穷多解! 自由变元个数为%d
", free_num);
148             for (i = 0; i < var; i++)
149             {
150                 if (free_x[i]) printf("x%d 是不确定的
", i + 1);
151                 else printf("x%d: %d
", i + 1, x[i]);
152             }
153         }
154         else
155         {
156             for (i = 0; i < var; i++)
157             {
158                 printf("x%d: %d
", i + 1, x[i]);
159             }
160         }
161         printf("
");
162     }
163     return 0;
164 }
165 
166 /* czyuan原创,转载请注明出处。*/
View Code

不过......我一般都是将这个精简过的代码改一下用(里面要用到求lcm最小公倍数和gcd最大共约数的函数)

 1 // 高斯消元法解方程组(Gauss-Jordan elimination).
 2 //未测试
 3 //equ表示方程数,var表示变量数
 4 //-2表示有浮点数解,但无整数解//-1表示无解//0表示唯一解//大于0表示无穷解,并返回自由变量的个数
 5 
 6 int a[N+5][N+5], x[N+5];        //a存储系数矩阵,x存储解
 7 
 8 int gauss(int equ, int var)
 9 {
10     int row = 0, col = 0;
11     while (row < equ && col < var){
12         int flag = row;
13         for (int i = row+1; i < equ; ++ i)
14             if (abs(a[i][col]) > abs(a[flag][col]))
15                 flag = i;
16         if (flag != row) for (int i = col; i <= var; ++ i)
17             swap (a[row][i], a[flag][i]);
18         if (!a[row][col]){ ++ col; continue;}
19 
20         for (int i = row+1; i < equ; ++ i){
21             if (a[i][col] == 0) continue;
22 
23             int tmp = lcm(abs(a[i][col]), abs(a[row][col]));
24             int ta = tmp / abs(a[i][col]);
25             int tb = tmp / abs(a[row][col]);
26             if (a[i][col] * a[row][col] < 0) tb = -tb;
27             for (int j = col; j <= var; ++ j){
28                 a[i][j] = (a[i][j] * ta - a[row][j]*tb) % 7;
29                 if (a[i][j] < 0) a[i][j] += 7;
30             }
31         }
32         ++ row; ++ col;
33     }
34     
35     for (int i = row; i < equ; ++ i)
36         if (a[i][var]) return -1;
37     if (row < var) return var - row;
38     for (int i = var-1; i >= 0; -- i){
39         int tmp = a[i][var];
40         for (int j = i+1; j < var; ++ j)
41             tmp -= a[i][j]*x[j];
42         if (tmp % a[i][i]) return -2;
43         x[i] = tmp / a[i][i];
44     }
45     return 0;
46 }
View Code

1、UVa 10828 Back to Kernighan-Ritchie

2、POJ 2947 Widget Factory

  1 /*
  2  * Author:  Plumrain
  3  * Created Time:  2013-10-07 16:41
  4  * File Name: math-POJ-2947.cpp
  5  */
  6 #include<iostream>
  7 #include<cstdio>
  8 #include<cstring>
  9 #include<cmath>
 10 #include<algorithm>
 11 
 12 using namespace std;
 13 
 14 #define CLR(x) memset(x, 0, sizeof(x))
 15 
 16 const int N = 300;
 17 
 18 int a[N+5][N+5], x[N+5];
 19 char ann[7] = {'O', 'U', 'E', 'H', 'R', 'A', 'U'};
 20 
 21 inline int gcd(int a, int b){
 22     return b == 0 ? a : gcd(b, a%b);
 23 }
 24 
 25 inline int lcm(int a, int b){
 26     return a * b / gcd(a, b);
 27 }
 28 
 29 int change(char a, char b)
 30 {
 31     for (int i = 0; i < 7; ++ i)
 32         if (i != 1 && i != 6 && ann[i] == b) return i;
 33     if (a == 'T') return 1;
 34     return 6;
 35 }
 36 
 37 void init(int equ, int var)
 38 {
 39     CLR (a);
 40     char s1[10], s2[10];
 41     int k;
 42     for (int i = 0; i < equ; ++ i){
 43         scanf ("%d%s%s", &k, s1, s2);
 44         a[i][var] = (change(s2[0], s2[1]) - change(s1[0], s1[1]) + 1) % 7;
 45         if (a[i][var] < 0) a[i][var] += 7;
 46         int tmp;
 47         while (k--){
 48             scanf ("%d", &tmp);
 49             ++ a[i][tmp-1];
 50             if (a[i][tmp-1] == 7) a[i][tmp-1] = 0;
 51         }
 52     }
 53 }
 54 
 55 int gauss(int equ, int var)
 56 {
 57     int row = 0, col = 0;
 58     while (row < equ && col < var){
 59         int flag = row;
 60         for (int i = row+1; i < equ; ++ i)
 61             if (abs(a[i][col]) > abs(a[flag][col]))
 62                 flag = i;
 63         if (flag != row) for (int i = col; i <= var; ++ i)
 64             swap (a[row][i], a[flag][i]);
 65         if (!a[row][col]){ ++ col; continue;}
 66 
 67         for (int i = row+1; i < equ; ++ i){
 68             if (a[i][col] == 0) continue;
 69 
 70             int tmp = lcm(abs(a[i][col]), abs(a[row][col]));
 71             int ta = tmp / abs(a[i][col]);
 72             int tb = tmp / abs(a[row][col]);
 73             if (a[i][col] * a[row][col] < 0) tb = -tb;
 74             for (int j = col; j <= var; ++ j){
 75                 a[i][j] = (a[i][j] * ta - a[row][j]*tb) % 7;
 76                 if (a[i][j] < 0) a[i][j] += 7;
 77             }
 78         }
 79         ++ row; ++ col;
 80     }
 81     
 82     for (int i = row; i < equ; ++ i)
 83         if (a[i][var]) return -1;
 84     if (row < var) return var - row;
 85     for (int i = var-1; i >= 0; -- i){
 86         int tmp = a[i][var];
 87         for (int j = i+1; j < var; ++ j)
 88             tmp = ((tmp - a[i][j]*x[j])%7 + 7) % 7;
 89         while (tmp % a[i][i] != 0) tmp += 7;
 90         x[i] = (tmp / a[i][i]) % 7;
 91     }
 92     return 0;
 93 }
 94 
 95 int main()
 96 {
 97     int equ, var;
 98     while (scanf ("%d%d", &var, &equ) != EOF && equ){
 99         init(equ, var);
100 
101         int ans = gauss(equ, var);
102         if (!ans){
103             for (int i = 0; i < var; ++ i)
104                 if (x[i] <= 2) x[i] += 7;
105             for (int i = 0; i < var-1; ++ i)
106                 printf ("%d ", x[i]);
107             printf ("%d
", x[var-1]);
108         }
109         else if (ans == -1)
110             printf ("Inconsistent data.
");
111         else
112             printf ("Multiple solutions.
");
113     }
114     return 0;
115 }
View Code

3、POJ 2065 SETI

  1 /*
  2  * Author:  Plumrain
  3  * Created Time:  2013-10-08 13:36
  4  * File Name: math-POJ-2065.cpp
  5  */
  6 #include<iostream>
  7 #include<cstdio>
  8 #include<cstring>
  9 #include<cmath>
 10 #include<algorithm>
 11 
 12 using namespace std;
 13 
 14 #define CLR(x) memset(x, 0, sizeof(x))
 15 
 16 const int N = 70;
 17 
 18 int a[N+5][N+5], mod, n, x[N+5];
 19 
 20 int gcd(int a, int b)
 21 {
 22     return b ? gcd(b, a%b) : a;
 23 }
 24 
 25 int lcm(int a, int b)
 26 {
 27     return a * b / gcd(a, b);
 28 }
 29 
 30 void init()
 31 {
 32     CLR (a); CLR (x);
 33     char s[N+5]; CLR (s);
 34     scanf ("%d%s", &mod, s);
 35     n = strlen(s);
 36     for (int i = 0; i < n; ++ i)
 37         if (s[i] != '*') a[i][n] = s[i] - 'a' + 1;
 38 
 39     for (int i = 0; i < n; ++ i){
 40         int tmp = 1, tim = 0;
 41         while (tim < n){
 42             a[i][tim++] = tmp;
 43             tmp = (tmp * (i+1)) % mod; 
 44         }
 45     }
 46 }
 47 
 48 int gauss()
 49 {
 50     int row = 0, col = 0;
 51     while (row < n && col < n){
 52         int flag = row;
 53         for (int i = row+1; i < n; ++ i)
 54             if (abs(a[i][col]) > abs(a[flag][col]))
 55                 flag = i;
 56         if (flag != row) for (int i = col; i <= n; ++ i)
 57             swap(a[flag][i], a[row][i]);
 58         if (!a[row][col]) {++ col; continue;}
 59         
 60         for (int i = row+1; i < n; ++ i){
 61             if (!a[i][col]) continue;
 62             
 63             int tmp = lcm(abs(a[i][col]), abs(a[row][col]));
 64             int ta = tmp / abs(a[i][col]);
 65             int tb = tmp / abs(a[row][col]);
 66             if (a[i][col] * a[row][col] < 0) tb = -tb;
 67             for (int j = col; j <= n; ++ j)
 68                 a[i][j] = ((a[i][j] * ta - a[row][j] * tb)%mod + mod) % mod;
 69         }
 70         ++ col; ++ row;
 71     }
 72 
 73     for (int i = n-1; i >= 0; -- i){
 74         int tmp = a[i][n];
 75         for (int j = i+1; j < n; ++ j)
 76             tmp = ((tmp - a[i][j] * x[j]) % mod + mod) % mod;
 77         if (a[i][i] < 0){a[i][i] = -a[i][i]; tmp = -tmp;}
 78         while (tmp % a[i][i]) tmp += mod;
 79         x[i] = (tmp / a[i][i]) % mod;
 80         if (x[i] < 0) x[i] += mod;
 81     }
 82     return 0;
 83 }
 84 
 85 int main()
 86 {
 87     int T;
 88     scanf ("%d", &T);
 89     while (T--){
 90         init(); 
 91 
 92         int ans = gauss();
 93 
 94         for (int i = 0; i < n; ++ i){
 95             if (i) printf (" ");
 96             printf ("%d", x[i]);
 97         }
 98         printf ("
");
 99     }
100     return 0;
101 }
View Code

以上三道题,感觉都是直接用模板解就好了,也没什么特别需要想的东西。

4、POJ 1166 The Clocks

这道题,感觉应该用搜索做的....看见网上有人说暴力可以过,有人说用高斯消元也行。但是个人感觉高斯消元应该不行啊,毕竟是要求最优解。而且,在用高斯消元的过程中碰到了一个问题,就是将当前要处理的行与别的行交换时,不能交换a[row][col]绝对值最大的行,也不能交换最后一个a[i][col]非零的行,而按照我下面的写法就恰好能过.....实在不明白原因,猜测这道题是卡数据过的- -。

 1 /*
 2  * Author:  Plumrain
 3  * Created Time:  2013-10-08 16:26
 4  * File Name: math-POJ-1166.cpp
 5  */
 6 #include<iostream>
 7 #include<cstdio>
 8 #include<cstring>
 9 
10 using namespace std;
11 
12 #define CLR(x) memset(x, 0, sizeof(x))
13 #define out(x) cout<<#x<<":"<<(x)<<endl
14 #define tst(a) cout<<#a<<endl
15 
16 int x[10], sum;
17 int a[10][10]={
18 };
19 
20 void gauss()
21 {
22     sum = 0;
23     int row = 0, col = 0;
24     while (row < 9 && col < 9){
25         int flag = row;
26         if (!a[row][col])
27             for (int i = row+1; i < 9; ++ i)
28                 if (a[i][col] && !a[flag][col]) flag = i;
29         if (flag != row) for (int i = col; i <= 9; ++ i)
30             swap (a[row][i], a[flag][i]);
31         if (!a[row][col]){++ col; continue;}
32 
33         for (int i = row+1; i < 9; ++ i){
34             int ta = a[i][col], tb = a[row][col];
35             for (int j = col; j <= 9; ++ j){
36                 a[i][j] = (a[i][j] * tb - a[row][j] * ta) % 4;
37                 if (a[i][j] < 0) a[i][j] += 4;
38             }
39         }
40         ++ col; ++ row;
41     }
42 
43     for (int i = 8; i >= 0; -- i){
44         int tmp = a[i][9];
45         for (int j = i+1; j < 9; ++ j)
46             tmp -= a[i][j] * x[j];
47         tmp %= 4; if (tmp < 0) tmp += 4;
48         for (x[i] = 0; x[i] < 4; ++ x[i])
49             if ((x[i]*a[i][i]%4 + 4) % 4 == tmp)
50                 break;
51         x[i] %= 4;
52         sum += x[i];
53     }
54 }
55 
56 int main()
57 {
58     CLR (a);
59     a[0][0]=a[1][0]=a[3][0]=a[4][0]=1;  
60     a[0][1]=a[1][1]=a[2][1]=1;  
61     a[1][2]=a[2][2]=a[4][2]=a[5][2]=1;  
62     a[0][3]=a[3][3]=a[6][3]=1;  
63     a[1][4]=a[3][4]=a[4][4]=a[5][4]=a[7][4]=1;  
64     a[2][5]=a[5][5]=a[8][5]=1;  
65     a[3][6]=a[4][6]=a[6][6]=a[7][6]=1;  
66     a[6][7]=a[7][7]=a[8][7]=1;  
67     a[4][8]=a[5][8]=a[7][8]=a[8][8]=1;  
68 
69     for (int i = 0; i < 9; ++ i){
70         scanf ("%d", &a[i][9]);
71         a[i][9] = (4 - a[i][9]) % 4;
72     }
73 
74     gauss();
75     
76     for (int i = 0; i < 9; ++ i)
77         while (x[i]){
78             printf ("%d", i+1);
79             -- x[i]; -- sum;
80             printf (sum>0 ? " " : "
");
81         }
82     return 0;
83 }
View Code

5、POJ 1222 EXTENDED LIGHTS OUT

题意:POJ的Discuss里面有人说可以证明此题只有唯一解。我不会证。此题可以用高斯消元的方法做,当然也可以用状态压缩暴力枚举第一行 + 递推的方式做。我比较偷懒,选择了后者。

 1 /*
 2  * Author:  Plumrain
 3  * Created Time:  2013-10-08 18:52
 4  * File Name: simulate-POJ-1222.cpp
 5  */
 6 #include<iostream>
 7 #include<cstdio>
 8 #include<cstring>
 9 #include<vector>
10 #include<utility>
11 
12 using namespace std;
13 
14 #define CLR(x) memset(x, 0, sizeof(x))
15 #define PB push_back
16 #define out(x) cout<<#x<<":"<<(x)<<endl
17 #define tst(a) cout<<#a<<endl
18 
19 int a[10][10], tmp[10][10], num[10][10];
20 vector<pair<int, int> > ans;
21 
22 void init()
23 {
24     ans.clear();
25     for (int i = 0; i < 5; ++ i)
26         for (int j = 0; j < 6; ++ j)
27             scanf ("%d", &a[i][j]);
28 
29     for (int i = 0; i < 5; ++ i)
30         for (int j = 0; j < 6; ++ j)
31             tmp[i][j] = a[i][j];
32 }
33 
34 void change(int x, int y)
35 {
36     ans.PB (make_pair(x, y));
37     a[x][y] = 1 - a[x][y];
38     if (x) a[x-1][y] = 1 - a[x-1][y];
39     if (x < 4) a[x+1][y] = 1 - a[x+1][y];
40     if (y) a[x][y-1] = 1 - a[x][y-1];
41     if (y < 5) a[x][y+1] = 1 - a[x][y+1];
42 }
43 
44 void all_init()
45 {
46     ans.clear();
47     for (int i = 0; i < 5; ++ i)
48         for (int j = 0; j < 6; ++ j)
49             a[i][j] = tmp[i][j];
50 }
51 
52 bool all_off(int x)
53 {
54     for (int i = 0; i < 6; ++ i)
55         if (a[x][i]) return false;
56     return true;
57 }
58 
59 void gao()
60 {
61     for (int i = 0; i < (1<<6); ++ i){
62         for (int j = 0; j < 6; ++ j)
63             if (i & (1<<j)) change(0, j);
64 
65         for (int j = 1; j < 5; ++ j)
66             for (int k = 0; k < 6; ++ k)
67                 if (a[j-1][k]) change(j, k);
68 
69         if (all_off(4)) return;
70         all_init();
71     }
72 }
73 
74 int main()
75 {
76     int T, test = 0;
77     scanf ("%d", &T);
78     while (T--){
79         printf ("PUZZLE #%d
", ++ test);
80 
81         init();
82         gao();
83 
84         CLR (num);
85         for (int i = 0; i < ans.size(); ++ i)
86             num[ans[i].first][ans[i].second] = 1;
87 
88         for (int i = 0; i < 5; ++ i){
89             for (int j = 0; j < 5; ++ j)
90                 printf ("%d ", num[i][j]);
91             printf ("%d
", num[i][5]);
92         }
93     }
94     return 0;
95 }
View Code

 

------------------------------------------------------------------
现在的你,在干什么呢?
你是不是还记得,你说你想成为岩哥那样的人。
原文地址:https://www.cnblogs.com/plumrain/p/gauss_elimination.html