高斯消元模板，整数（数学）

Posted 2020-06-14

线性代数，原理不讲了。。。

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