跳转至

高斯消元⚓︎

高斯消元(\text{Gauss Elimination})是一种用于解线性方程组的算法。它通过对增广矩阵进行行变换,将其转化为上三角矩阵,然后通过回代求解变量的值。

高斯消元法可以解线性方程组、异或方程组、同余方程组等问题。
高斯消元法的时间复杂度为 O(n^3),其中 n 是方程组中变量的数量。

线性方程组⚓︎

线性方程组可以表示为增广矩阵的形式,通过高斯消元法求解。

线性方程组
C++
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <utility>
#include <vector>
using namespace std;

const double eps = 1e-9;

void GaussElimination(vector<vector<double>> &matrix) {
  int n = matrix.size();
  // 消元
  for (int i = 0; i < n; ++i) {
    // 寻找主元, 为了区分多解和无解的情况, 这里需要从第0行开始找
    int pivot = i;
    for (int j = 0; j < n; ++j) {
      // 已经消元的行不再考虑
      if (j < i && abs(matrix[j][j]) >= eps) { continue; }
      if (abs(matrix[j][i]) > abs(matrix[pivot][i])) { pivot = j; }
    }
    swap(matrix[i], matrix[pivot]);

    if (abs(matrix[i][i]) < eps) { continue; }  // 主元为0, 无法消元
    // 将主元化为1
    double major = matrix[i][i];
    for (int j = i; j < n + 1; ++j) { matrix[i][j] /= major; }
    // 消元
    for (int j = 0; j < n; ++j) {
      if (j != i) {
        double factor = matrix[j][i];
        for (int k = i; k < n + 1; ++k) { matrix[j][k] -= factor * matrix[i][k]; }
      }
    }
  }
}

void BackSubstitution(const vector<vector<double>> &matrix) {
  int n = matrix.size();
  // 回代
  vector<double> result(n, 0);
  bool multiple_solutions = false;
  for (int i = 0; i < n; ++i) {
    if (abs(matrix[i][i]) < eps) {
      // 无解: 该行主元为0, 但常数项不为0
      if (abs(matrix[i][n]) >= eps) {
        cout << -1 << '\n';
        return;
      }
      // 多解, 该行主元为0, 且常数项也为0 (自由元). 注意后续可能出现无解的情况
      multiple_solutions = true;
    }
    result[i] = matrix[i][n];
  }
  if (multiple_solutions) {  // 存在无穷多解
    cout << 0 << '\n';
    return;
  }
  if (!result.empty()) {
    int count = 1;
    for (double v : result) {
      cout << "x" << count++ << "=" << fixed << setprecision(2) << v << "\n";
    }
  }
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  cout.tie(nullptr);
  int n;
  cin >> n;
  // 输入增广矩阵: n*(n+1)
  vector<vector<double>> matrix(n, vector<double>(n + 1));
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) { cin >> matrix[i][j]; }
    cin >> matrix[i][n];
  }
  GaussElimination(matrix);
  BackSubstitution(matrix);
  return 0;
}

异或方程组⚓︎

异或(\text{XOR})方程组是一类特殊的线性方程组,其中变量和常数项都在 GF(2)(即模 2 的有限域)上进行运算。
异或方程组可以通过高斯消元法在 GF(2) 上求解。

Lights G

n 个开关和 n 个灯泡,每个开关控制一个灯泡的开关状态(开或关)。
每按下一个开关,除了该开关控制的灯泡状态会改变外,与该灯泡相邻的灯泡状态也会改变。
现在给定初始状态和目标状态,问最少按下多少次开关可以达到目标状态。

C++
#include <functional>
#include <iostream>
#include <vector>
using namespace std;

bool GaussXor(vector<vector<int>> &matrix) {
  int n = matrix.size();
  // 消元
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      if (j < i && matrix[j][j] == 1) { continue; }
      if (matrix[j][i] == 1) {
        swap(matrix[i], matrix[j]);
        break;
      }
    }
    if (matrix[i][i] == 0) { continue; }
    for (int j = 0; j < n; ++j) {
      if (j != i && matrix[j][i] == 1) {
        for (int k = i; k < n + 1; ++k) { matrix[j][k] ^= matrix[i][k]; }
      }
    }
  }

  for (int i = 0; i < n; ++i) {
    if (matrix[i][i] == 0) { return false; }
  }
  return true;
}

bool BackSubstitution(const vector<vector<int>> &matrix) {
  int n = matrix.size();
  vector<int> result(n, 0);
  bool multiple_solutions = false;
  for (int i = 0; i < n; ++i) {
    if (matrix[i][i] == 0) {
      // 无解: 该行主元为0, 但常数项不为0
      if (matrix[i][n] == 1) { return false; }
      // 多解, 该行主元为0, 且常数项也为0 (自由元). 注意后续可能出现无解的情况
      multiple_solutions = true;
    }
    result[i] = matrix[i][n];
  }
  // 多解
  if (multiple_solutions) { return true; }
  // 唯一解
  return true;
}

int main() {
  int n, m;
  cin >> n >> m;
  // 构造增广矩阵: n*(n+1)
  vector<vector<int>> matrix(n, vector<int>(n + 1));
  for (int i = 0; i < n; ++i) {
    matrix[i][i] = 1;
    matrix[i][n] = 1;
  }
  for (int i = 0; i < m; ++i) {
    int a, b;
    cin >> a >> b;
    --a, --b;
    matrix[a][b] = 1;
    matrix[b][a] = 1;
  }
  int ans           = 0;
  bool has_solution = GaussXor(matrix);
  if (has_solution) {
    for (int i = 0; i < n; ++i) { ans += matrix[i][n]; }
  } else {
    ans = n;
    vector<int> ops(n, 0);
    std::function<void(int, int)> dfs = [&](int i, int res) {
      if (res >= ans) { return; }
      if (i == -1) {
        ans = min(ans, res);
        return;
      }
      if (matrix[i][i] == 0) {  // 自由元
        ops[i] = 0;
        dfs(i - 1, res);
        ops[i] = 1;
        dfs(i - 1, res + 1);
      } else {  // 主元
        int val = matrix[i][n];
        for (int j = i + 1; j < n; ++j) {
          if (matrix[i][j] == 1) { val ^= ops[j]; }
        }
        dfs(i - 1, res + val);
      }
    };
    dfs(n - 1, 0);
  }
  cout << ans << "\n";
  return 0;
}

同余方程组⚓︎

同余方程组是一类特殊的线性方程组,其中变量和常数项都在模 m 的整数环上进行运算。
同余方程组可以通过扩展的欧几里得算法和中国剩余定理求解。
高斯消元法也可以用于求解同余方程组,但需要对行变换进行适当的调整。

同余方程组
C++
#include <numeric>
#include <vector>
using namespace std;

void GaussModular(vector<vector<int>> &matrix, int64_t mod) {
  int n = matrix.size();
  // 增广矩阵: n*(n+1), 处理系数矩阵, 使其非负
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j <= n; ++j) { matrix[i][j] = (matrix[i][j] % mod + mod) % mod; }
  }

  // 消元
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      // 已经消元的行不再考虑
      if (j < i && matrix[j][j] != 0) { continue; }
      if (matrix[j][i] != 0) {
        swap(matrix[i], matrix[j]);
        break;
      }
    }

    if (matrix[i][i] == 0) { continue; }  // 主元为0, 无法消元
    for (int j = 0; j < n; ++j) {
      if (j != i && matrix[j][i] != 0) {
        int factor   = gcd(matrix[j][i], matrix[i][i]);
        int factor_i = matrix[i][i] / factor;
        int factor_j = matrix[j][i] / factor;
        if (j < i && matrix[j][j] != 0) {
          // 如果j行有主元,那么从j列到i-1列的所有系数 * factor_j
          // 正确更新主元和自由元之间的关系
          for (int k = j; k < i; ++k) { matrix[j][k] = (matrix[j][k] * factor_i) % mod; }
        }
        // 消元
        for (int k = i; k < n + 1; ++k) {
          matrix[j][k] = ((matrix[j][k] * factor_i - matrix[i][k] * factor_j) % mod + mod) % mod;
        }
      }
    }
  }
  // 逆元表, 只适用于质数模
  vector<int64_t> inverse(mod, 0);
  inverse[1] = 1;
  for (int i = 2; i < mod; ++i) { inverse[i] = (mod - mod / i) * inverse[mod % i] % mod; }
  // 收尾, 设置主元为1
  for (int i = 0; i < n; ++i) {
    if (matrix[i][i] == 0) { continue; }
    bool free_variable = false;
    for (int j = i + 1; j < n; ++j) {
      if (matrix[i][j] != 0) {
        free_variable = true;
        break;
      }
    }
    // 如果该行有自由元,不能直接归一化主元
    if (free_variable) { continue; }
    int64_t inv  = inverse[matrix[i][i]];
    matrix[i][i] = 1;
    matrix[i][n] = (matrix[i][n] * inv) % mod;
  }
}

// 回代
void BackSubstitution(const vector<vector<double>> &matrix) {
  int n = matrix.size();
  vector<int> result(n, 0);
  bool multiple_solutions = false;
  for (int i = 0; i < n; ++i) {
    if (matrix[i][i] == 0) {
      // 无解: 该行主元为0, 但常数项不为0
      if (matrix[i][n] != 0) { return; }
      // 多解, 该行主元为0, 且常数项也为0 (自由元). 注意后续可能出现无解的情况
      multiple_solutions = true;
    }
    result[i] = matrix[i][n];
  }
  // 多解
  if (multiple_solutions) { return; }
  return;
}
Gambler Bo

nm 列的格子,每个格子有一个初始值。
每次操作可以选择一个格子,令该格子的值加 2 和它上下左右四个格子的值加 1(如果存在)。
目标是通过若干次操作使得所有格子的值都能被 3 整除(即模 3 等于 0)。

数据范围

T \leq 10, 1 \leq N,M \leq 30
如果每次消元从 0 行开始找主元,可能会超时。

C++
#include <array>
#include <iostream>
#include <utility>
#include <vector>
using namespace std;

const int MOD = 3;

// 求模 3 下逆元
inline int inv3(int x) {
  if (x == 1) { return 1; }
  if (x == 2) { return 2; }
  return 0;  // 0 无逆元
}

// 模 3 高斯消元
// matrix: n x (n+1) 增广矩阵
// x: 返回解
void GaussModular(vector<vector<int>> &matrix, vector<int> &x, int n) {
  for (int i = 0; i < n; ++i) {
    // 找主元
    int pivot = -1;
    for (int j = i; j < n; ++j) {
      if (matrix[j][i] != 0) {
        pivot = j;
        break;
      }
    }
    if (pivot == -1) { continue; }  // 主元为0, 无法消元
    swap(matrix[i], matrix[pivot]);

    // 主元归一化
    int inv = inv3(matrix[i][i]);
    for (int k = i; k <= n; ++k) { matrix[i][k] = (matrix[i][k] * inv) % MOD; }

    // 消元整列
    for (int j = 0; j < n; ++j) {
      if (j == i || matrix[j][i] == 0) { continue; }
      int factor = matrix[j][i];
      for (int k = i; k <= n; ++k) {
        matrix[j][k] = (matrix[j][k] - factor * matrix[i][k] % MOD + MOD) % MOD;
      }
    }
  }

  // 回代
  for (int i = 0; i < n; ++i) { x[i] = matrix[i][n]; }
}

// 回代输出结果
void OutputSolution(const vector<int> &x, int n, int m) {
  int total = 0;
  for (int i = 0; i < n * m; ++i) { total += x[i]; }
  cout << total << "\n";
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      int idx = i * m + j;
      for (int k = 0; k < x[idx]; ++k) { cout << i + 1 << " " << j + 1 << "\n"; }
    }
  }
}

// 构建增广矩阵
void BuildMatrix(vector<vector<int>> &matrix, int n, int m, const vector<int> &vals) {
  array<int, 4> dir_x = {0, 0, 1, -1};
  array<int, 4> dir_y = {1, -1, 0, 0};
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      int idx          = i * m + j;
      matrix[idx][idx] = 2;  // 对自己影响为 2
      // 对四个方向影响为 1
      for (int d = 0; d < 4; ++d) {
        int ni = i + dir_x[d];
        int nj = j + dir_y[d];
        if (ni < 0 || ni >= n || nj < 0 || nj >= m) { continue; }
        matrix[idx][ni * m + nj] = 1;
      }
      matrix[idx][n * m] = (MOD - vals[idx]) % MOD;
    }
  }
}

void solve() {
  int n, m;
  cin >> n >> m;
  vector<int> vals(n * m);
  for (int i = 0; i < n * m; ++i) { cin >> vals[i]; }

  vector<vector<int>> matrix(n * m, vector<int>(n * m + 1, 0));
  vector<int> x(n * m, 0);
  BuildMatrix(matrix, n, m, vals);
  GaussModular(matrix, x, n * m);
  OutputSolution(x, n, m);
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  cout.tie(nullptr);
  int t;
  cin >> t;
  while ((t--) != 0) { solve(); }
  return 0;
}

评论