zoukankan      html  css  js  c++  java
  • 高斯消元问题

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

    这里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

     

    ------------------------------------------------------------------
    现在的你,在干什么呢?
    你是不是还记得,你说你想成为岩哥那样的人。
  • 相关阅读:
    NTP服务器搭建
    Linux安装MongoDB 4.4.2
    CentOS安装Zookeeper 3.6.2
    CentOS安装Redis 6.0.9
    MacBook Home End
    SLES Install
    cucumber soapui test web services
    S/4 HANA Solution Manager
    Linux下创建新用户
    su with hyphen and without
  • 原文地址:https://www.cnblogs.com/plumrain/p/gauss_elimination.html
Copyright © 2011-2022 走看看