高斯消元小结

Author Avatar
TaoSama Aug 05, 2016

Ⅰ. 高斯消元认知

  • $这东西英文名叫Gaussian Elimination$
  • $高斯消元对矩阵进行操作，对满足一定运算规律的方程利用矩阵初等变换进行消元$
  • $主要用来求解线性方程组、矩阵的秩、以及可逆方阵的逆矩阵$
  • $还可以用来线性空间的基向量，比如线性基（线性无关的极大子集）$
  • $PS:关于线性相关，具体可以读2014集训队论文《浅谈线性相关》$

Ⅱ. 高斯消元模版

浮点数高斯消元

double a[N][N], ans[N];
bool isFreeX[N];

void getAns(int n, int m, int r) {
//    memset(ans, 0, sizeof ans);  //not necessary
    for(int i = r - 1; ~i; --i) {
        for(int j = 0; j < m; ++j) {
            if(!sgn(a[i][j])) continue;
            ans[j] = a[i][m];
            for(int k = j + 1; k < m; ++k)
                ans[j] -= a[i][k] * ans[k];
            ans[j] = ans[j] / a[i][j];
            break;
        }
    }
}

int gauss(int n, int m) {
    for(int i = 0; i < m; ++i) isFreeX[i] = false;
    int r = 0, c = 0;
    for(; r < n && c < m; ++r, ++c) {
        int maxR = r;       //row transform
        for(int i = r + 1; i < n; ++i)
            if(abs(a[i][c]) > abs(a[maxR][c])) maxR = i;
        if(maxR != r) swap(a[maxR], a[r]);

        if(!sgn(a[r][c])) { --r; isFreeX[c] = true; continue;}

        //eliminate coefficient
        for(int i = r + 1; i < n; ++i) {
            if(sgn(a[i][c])) {
                double delta = a[i][c] / a[r][c];
                for(int j = c; j <= m; ++j)
                    a[i][j] -= delta * a[r][j];
            }
        }
    }
    for(int i = r; i < n; i++) if(sgn(a[i][m])) return -1;

    getAns(n, m, r);

    //at last, r is rank, m - r is the number of freeX
    return r;
}

模意义下高斯消元

int a[N][N], ans[N];
bool isFreeX[N];

inline int inv(int x) {return 1;}

int getAns(int n, int m, int r) {
    int ret = 0;
    for(int i = r - 1; ~i; --i) {
        for(int j = 0; j < m; ++j) {
            if(!a[i][j]) continue;
            ans[j] = a[i][m];
            for(int k = j + 1; k < m; ++k) {
                ans[j] -= a[i][k] * ans[k];
                ans[j] %= MOD;
                if(ans[j] < 0) ans[j] += MOD;
            }
            ans[j] = ans[j] * inv(a[i][j]) % MOD;
            break;
        }
    }
    for(int i = 0; i < m; ++i) ret += ans[i];
    return ret;
}

int gauss(int n, int m) {
    for(int i = 0; i < m; ++i) isFreeX[i] = false;
    int r = 0, c = 0;
    for(; r < n && c < m; ++r, ++c) {
        int maxR = r;       //row transform
        for(int i = r + 1; i < n; ++i)
            if(abs(a[i][c]) > abs(a[maxR][c])) maxR = i;
        if(maxR != r) swap(a[maxR], a[r]);

        if(!a[r][c]) { --r; isFreeX[c] = true; continue;}

        //eliminate coefficient
        for(int i = r + 1; i < n; ++i) {
            if(a[i][c]) {
                int delta = a[i][c] * inv(a[r][c]);
                for(int j = c; j <= m; ++j) {
                    a[i][j] -= delta * a[r][j];
                    a[i][j] %= MOD;
                    if(a[i][j] < 0) a[i][j] += MOD;
                }
            }
        }
    }
    for(int i = r; i < n; i++) if(a[i][m]) return -1;

    //at last, r is rank, m - r is the number of freeX
    return r;
}

int getMinAns(int n, int m, int r) {
    int ret = INF, freeX = m - r;
    for(int s = 0; s < 1 << freeX; ++s) {
        if(__builtin_popcount(s) >= ret) continue;
        int cnt = 0;
        for(int j = 0; j < m; ++j) {
            if(isFreeX[j]) {
                ans[j] = s >> cnt & 1;
                ++cnt;
            }
        }
        ret = min(ret, getAns(n, m, r));
    }
    return ret;
}

XOR高斯消元

int a[N][N], ans[N];

int xorGauss() {
    int r = 0, c = 0;
    for(; r < n && c < n; ++r, ++c) {
        int p = r;
        for(; p < n; ++p) if(a[p][c]) break;
        if(p == n) {--r; continue;}
        swap(a[p], a[r]);

        for(int i = 0; i < n; ++i) {
            if(i != r && a[i][c]) {
                for(int j = c; j <= n; ++j)
                    a[i][j] ^= a[r][j];
            }
        }
    }
    for(int i = r; i < n; ++i) if(a[i][n]) return -1;

    //memset(ans, 0, sizeof ans);  //not necessary
    for(int i = r - 1; ~i; --i)
        for(int j = 0; j < n; ++j)
            if(a[i][j]) {ans[j] = a[i][n]; break;}
    return r;
}

Ⅲ. 题目讲解

HDU 4200 Bad Wiring

分析:
$开关问题，POJ的几个数据都比较弱，所以就选了这个$
$建立模2意义下的方程组，然后高斯消元，要求最小解，二进制枚举自由变元回代即可$
代码：

//
//  Created by TaoSama on 2016-07-28
//  Copyright (c) 2016 TaoSama. All rights reserved.
//
#pragma comment(linker, "/STACK:102400000,102400000")
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <string>
#include <set>
#include <vector>

using namespace std;
#define pr(x) cout << #x << " = " << x << "  "
#define prln(x) cout << #x << " = " << x << endl
const int N = 100 + 10, INF = 0x3f3f3f3f, MOD = 2;

int n, d;
int a[N][N], ans[N];
bool isFreeX[N];

inline int inv(int x) {return 1;}

int getAns(int n, int m, int r) {
    int ret = 0;
    for(int i = r - 1; ~i; --i) {
        for(int j = 0; j < m; ++j) {
            if(!a[i][j]) continue;
            ans[j] = a[i][m];
            for(int k = j + 1; k < m; ++k) {
                ans[j] -= a[i][k] * ans[k];
                ans[j] %= MOD;
                if(ans[j] < 0) ans[j] += MOD;
            }
            ans[j] = ans[j] * inv(a[i][j]) % MOD;
            break;
        }
    }
    for(int i = 0; i < m; ++i) ret += ans[i];
    return ret;
}

int gauss(int n, int m) {
    for(int i = 0; i < m; ++i) isFreeX[i] = false;
    int r = 0, c = 0;
    for(; r < n && c < m; ++r, ++c) {
        int maxR = r;       //row transform
        for(int i = r + 1; i < n; ++i)
            if(abs(a[i][c]) > abs(a[maxR][c])) maxR = i;
        if(maxR != r) swap(a[maxR], a[r]);

        if(!a[r][c]) { --r; isFreeX[c] = true; continue;}

        //eliminate coefficient
        for(int i = r + 1; i < n; ++i) {
            if(a[i][c]) {
                int delta = a[i][c] * inv(a[r][c]);
                for(int j = c; j <= m; ++j) {
                    a[i][j] -= delta * a[r][j];
                    a[i][j] %= MOD;
                    if(a[i][j] < 0) a[i][j] += MOD;
                }
            }
        }
    }
    for(int i = r; i < n; i++) if(a[i][m]) return -1;

    //at last, r is rank, m - r is the number of freeX
    return r;
}

int getMinAns(int n, int m, int r) {
    int ret = INF, freeX = m - r;
    for(int s = 0; s < 1 << freeX; ++s) {
        if(__builtin_popcount(s) >= ret) continue;
        int cnt = 0;
        for(int j = 0; j < m; ++j) {
            if(isFreeX[j]) {
                ans[j] = s >> cnt & 1;
                ++cnt;
            }
        }
        ret = min(ret, getAns(n, m, r));
    }
    return ret;
}

int main() {
#ifdef LOCAL
    freopen("C:\\Users\\TaoSama\\Desktop\\in.txt", "r", stdin);
//  freopen("C:\\Users\\TaoSama\\Desktop\\out.txt","w",stdout);
#endif
    ios_base::sync_with_stdio(0);
    clock_t _ = clock();

    int t; scanf("%d", &t);
    while(t--) {
        scanf("%d%d", &n, &d);
        memset(a, 0, sizeof a);
        for(int i = 0; i < n; ++i) {
            int v; scanf("%d", &v);
            a[i][n] = -v + MOD;
            for(int j = max(0, i - d); j <= min(n - 1, i + d); ++j)
                a[i][j] = 1;
        }
        int r = gauss(n, n);
        if(~r) printf("%d\n", getMinAns(n, n, r));
        else puts("impossible");
    }

#ifdef LOCAL
    printf("\nTime cost: %.2fs\n", 1.0 * (clock() - _) / CLOCKS_PER_SEC);
#endif
    return 0;
}

SGU 275 To xor or not to xor

分析：
$求N个数的子集的最大异或和，高斯消元后求出的是极大线性无关组即线性基$
$线性基异或起来显然就是最大异或和了$
代码：

//
//  Created by TaoSama on 2016-07-28
//  Copyright (c) 2016 TaoSama. All rights reserved.
//
#pragma comment(linker, "/STACK:102400000,102400000")
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <string>
#include <set>
#include <vector>

using namespace std;
#define pr(x) cout << #x << " = " << x << "  "
#define prln(x) cout << #x << " = " << x << endl
const int N = 1e5 + 10, INF = 0x3f3f3f3f, MOD = 1e9 + 7;

typedef long long LL;
int n;
LL a[N];

int xorGauss() {
    int r = 0, c = 63;
    for(; r < n && ~c; ++r, --c) {
        int p = r;
        for(; p < n; ++p) if(a[p] >> c & 1) break;
        if(p == n) {--r; continue;}
        swap(a[p], a[r]);

        for(int i = 0; i < n; ++i)
            if(i != r && a[i] >> c & 1) a[i] ^= a[r];
    }
    return r;
}

int main() {
#ifdef LOCAL
    freopen("C:\\Users\\TaoSama\\Desktop\\in.txt", "r", stdin);
//  freopen("C:\\Users\\TaoSama\\Desktop\\out.txt","w",stdout);
#endif
    ios_base::sync_with_stdio(0);
    clock_t _ = clock();

    scanf("%d", &n);
    for(int i = 0; i < n; ++i) scanf("%lld", a + i);

    int r = xorGauss();
    LL ans = 0;
    for(int i = 0; i < r; ++i) ans ^= a[i];
    printf("%lld\n", ans);


#ifdef LOCAL
    printf("\nTime cost: %.2fs\n", 1.0 * (clock() - _) / CLOCKS_PER_SEC);
#endif
    return 0;
}

HDU 2262 Where is the canteen

分析：
$高斯消元求解期望dp$
$f[x][y]:=(x, y)到终点的期望次数，显然终点的期望是0，即f[ex][ey]=0$
$期望倒推，f[x][y]=(\sum p_i\times f[nx][ny]) + 1$
$bfs求出每个点的可达性，之后对于每个点建立期望方程，然后高斯消元求解即可$
代码：

//
//  Created by TaoSama on 2016-07-28
//  Copyright (c) 2016 TaoSama. All rights reserved.
//
#pragma comment(linker, "/STACK:102400000,102400000")
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <string>
#include <set>
#include <vector>

using namespace std;
#define pr(x) cout << #x << " = " << x << "  "
#define prln(x) cout << #x << " = " << x << endl
const int N = 15 * 15 + 10, INF = 0x3f3f3f3f, MOD = 1e9 + 7;
const double EPS = 1e-8;

int sgn(double x) {
    return x < -EPS ? -1 : x > EPS;
}

int n, m;
char s[N][N];
double a[N][N], ans[N];
bool isFreeX[N];

void getAns(int n, int m, int r) {
//    memset(ans, 0, sizeof ans);  //not necessary
    for(int i = r - 1; ~i; --i) {
        for(int j = 0; j < m; ++j) {
            if(!sgn(a[i][j])) continue;
            ans[j] = a[i][m];
            for(int k = j + 1; k < m; ++k)
                ans[j] -= a[i][k] * ans[k];
            ans[j] = ans[j] / a[i][j];
            break;
        }
    }
}

int gauss(int n, int m) {
    for(int i = 0; i < m; ++i) isFreeX[i] = false;
    int r = 0, c = 0;
    for(; r < n && c < m; ++r, ++c) {
        int maxR = r;       //row transform
        for(int i = r + 1; i < n; ++i)
            if(abs(a[i][c]) > abs(a[maxR][c])) maxR = i;
        if(maxR != r) swap(a[maxR], a[r]);

        if(!sgn(a[r][c])) { --r; isFreeX[c] = true; continue;}

        //eliminate coefficient
        for(int i = r + 1; i < n; ++i) {
            if(sgn(a[i][c])) {
                double delta = a[i][c] / a[r][c];
                for(int j = c; j <= m; ++j)
                    a[i][j] -= delta * a[r][j];
            }
        }
    }
    for(int i = r; i < n; i++) if(sgn(a[i][m])) return -1;

    getAns(n, m, r);

    //at last, r is rank, m - r is the number of freeX
    return r;
}

inline int ID(int x, int y) {
    return x * m + y;
}

bool can[N][N];
const int d[][2] = { -1, 0, 0, -1, 1, 0, 0, 1};
bool bfs(int sx, int sy) {
    memset(can, 0, sizeof can);
    queue<pair<int, int> > q;
    q.push({sx, sy});
    can[sx][sy] = true;

    bool ok = false;
    while(q.size()) {
        auto u = q.front(); q.pop();
        int x, y; tie(x, y) = u;
        if(s[x][y] == '$') ok = true;
        for(int i = 0; i < 4; ++i) {
            int nx = x + d[i][0], ny = y + d[i][1];
            if(nx < 0 || nx >= n || ny < 0 || ny >= m) continue;
            if(s[nx][ny] == '#' || can[nx][ny]) continue;
            can[nx][ny] = true;
            q.push({nx, ny});
        }
    }
    return ok;
}

int main() {
#ifdef LOCAL
    freopen("C:\\Users\\TaoSama\\Desktop\\in.txt", "r", stdin);
//  freopen("C:\\Users\\TaoSama\\Desktop\\out.txt","w",stdout);
#endif
    ios_base::sync_with_stdio(0);
    clock_t _ = clock();

    while(scanf("%d%d", &n, &m) == 2) {
        int sx, sy; sx = sy = -1;
        for(int i = 0; i < n; ++i) {
            scanf("%s", s[i]);
            if(~sx) continue;
            for(int j = 0; j < m; ++j)
                if(s[i][j] == '@') sx = i, sy = j;
        }
        if(!bfs(sx, sy)) {puts("-1"); continue;}

        //E_i = sum(E_j)/cnt + 1 => sum(E_j) - cnt*E_i = -cnt
        memset(a, 0, sizeof a);
        for(int i = 0; i < n; ++i) {
            for(int j = 0; j < m; ++j) {
                if(s[i][j] == '#') continue; //obstacle
                if(s[i][j] == '$') { //destination
                    a[ID(i, j)][ID(i, j)] = 1;
                    a[ID(i, j)][n * m] = 0; //E_exit = 0;
                    continue;
                }

                int cnt = 0;
                for(int k = 0; k < 4; ++k) {
                    int x = i + d[k][0], y = j + d[k][1];
                    if(x < 0 || x >= n || y < 0 || y >= m || !can[x][y]) continue;
                    ++cnt;
                    a[ID(i, j)][ID(x, y)] = 1;
                }
                a[ID(i, j)][ID(i, j)] = -cnt;
                a[ID(i, j)][n * m] = -cnt;
            }
        }

//        for(int i = 0; i < n * m; ++i)
//          for(int j = 0; j <= n * m; ++j)
//              printf("%.3f%c", a[i][j], " \n"[j == n * m]);

        gauss(n * m, n * m);
        double result = ans[ID(sx, sy)];
        printf("%.6f\n", result);
    }

#ifdef LOCAL
    printf("\nTime cost: %.2fs\n", 1.0 * (clock() - _) / CLOCKS_PER_SEC);
#endif
    return 0;
}