高斯消元解线性方程组

输入一个包含 n 个方程 n 个未知数的线性方程组，方程组中的系数为实数。求解这个方程组。

题目链接：AcWing 883。

通过初等行变换转化为上三角矩阵，并回代，得到方程组的解。

对于方程组：

{a11x1+⋯+a1nxn=b1a21x1+⋯+a2nxn=b2⋮an1x1+⋯+annxn=bn\left \{ \begin{aligned} &a_{11}x_1+\cdots+a_{1n}x_n=b_1\\ &a_{21}x_1+\cdots+a_{2n}x_n=b_2\\ &\vdots\\ &a_{n1}x_1+\cdots+a_{nn}x_n=b_n \end{aligned} \right .

写成矩阵的形式就是：

Ax=bAx=b

系数矩阵 AA 的增广矩阵 [A∣b][A|b] 是系数矩阵加了一列右边的常数：

[a11⋯b1⋮⋱⋮an1⋯bn]\begin{bmatrix} a_{11}&\cdots&b_1\\ \vdots&\ddots&\vdots\\ a_{n1}&\cdots&b_n \end{bmatrix}

消元的目标就是将其变为一个上三角矩阵，然后再把它由下到上变换为一个类似单位矩阵的形式：

[1t11t2⋱⋮1tn−11tn]\begin{bmatrix} 1&&&&&t_1\\ &1&&&&t_2\\ &&\ddots&&&\vdots\\ &&&1&&t_{n-1}\\ &&&&1&t_n \end{bmatrix}

消元后的矩阵是 [I∣A−1b][I|A^{-1}b]，由于 x=A−1bx=A^{-1}b，最后一列的 t1,⋯ ,tnt_1,\cdots,t_n 就是方程组的解。

当最终的矩阵的对角线上出现 0，若对应解为 0，那么方程组有无穷多个解；如果对应解不为 0，方程组无解。

算法的步骤如下：

  1. 为了保持精度，先找到当前主元绝对值最大的一行，与当前行交换。
  2. 用主元除该行元素。
  3. 将主元下面同一列的元素全部通过初等行列变换为 0。
  4. 最后消掉上三角中的元素。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <cstdio>
#include <cmath>
using namespace std;

const int N = 110;
const double eps = 1e-6;
double a[N][N];
int n;

// 0: 有解, 1: 无解, 2: 无穷多解
int gauss() {
    // p: pivot, r: row
    int p = 0, r = 0;
    for (; p < n; p++) {
        // 寻找当前主元绝对值最大的一行
        int t = r;
        for (int i = r; i < n; i++)
            if (fabs(a[i][p]) > fabs(a[t][p]))
                t = i;
        
        // 主元为 0 跳过
        if (fabs(a[t][p]) < eps) continue;
        
        // 交换绝对值最大的主元所在行 与当前行
        for (int j = p; j <= n; j++)
            swap(a[r][j], a[t][j]);
        
        // 用主元除当前行 这里要逆序否则就要多加一个临时变量了
        for (int j = n; j >= p; j--)
            a[r][j] /= a[r][p];
        
        // 用主元所在行消去其所在列
        for (int i = r+1; i < n; i++)
            // if (fabs(a[i][p]) > eps)
            for (int j = n; j >= p; j--)
                a[i][j] -= a[i][p] * a[r][j];
        
        r++;
    }
    // 不够 n 个方程
    if (r <	n) {
        // 主元为 0 的行都被换下去了 r~n-1
        for (int i = r; i < n; i++)
            if (fabs(a[i][n]) > eps)
                return 2;
        return 1;
    }
    
    for (int i = n-1; i >= 0; i--)
        for (int j = i+1; j < n; j++)
            // a[j][n] 是 xj 的解, a[i][j] 是 xj 的系数
            a[i][n] -= a[i][j] * a[j][n];
    return 0;
}

int main() {
    scanf("%d", &n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j <= n; j++)
            scanf("%lf", &a[i][j]);
    int t = gauss();
    if (t == 1) puts("Infinite group solutions");
    else if (t == 2) puts("No solution");
    else {
        for (int i = 0; i < n; i++)
            printf("%.2lf\n", a[i][n]);
    }
    return 0;
}

球形空间产生器

在 nn 维空间的球面上给定 n+1n+1 个点的坐标，求球心坐标。数据保证唯一解。

题目链接：AcWing 207，P4035。

本题也可以用模拟退火做，但这里选择用高斯消元做，能保证求出精确解。

设球心坐标 (x1,x2,⋯ ,xn)(x_1,x_2,\cdots,x_n) 列出这 n+1n+1 个点的方程：

{(x1−a11)2+⋯+(xn−a1n)2=R2(x1−a21)2+⋯+(xn−a2n)2=R2⋮(x1−an+1,1)2+⋯+(xn−an+1,n)2=R2\left\{ \begin{aligned} &(x_1-a_{11})^2+\cdots+(x_n-a_{1n})^2=R^2\\ &(x_1-a_{21})^2+\cdots+(x_n-a_{2n})^2=R^2\\ &\vdots\\ &(x_1-a_{n+1,1})^2+\cdots+(x_n-a_{n+1,n})^2=R^2 \end{aligned} \right.

但这里是一个二次方程组，高斯消元能解决的是一次方程组，所以这里用每一个式子减去上一个式子，得到 nn 个线性方程。用 2 式减 1 式为例：

∑i=1n2(a2i−a1i)xi=∑i=1n(a2i2−a1i2)\sum_{i=1}^n2(a_{2i}-a_{1i})x_i = \sum_{i=1}^n(a_{2i}^2-a_{1i}^2)\\

这样就能构造出 nn 个线性方程组了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const int N = 15;
int n;
double a[N][N];

void gauss() {
    for (int r = 0; r < n; r++) {
        int t = r;
        for (int i = r; i < n; i++)
            if (fabs(a[i][r]) > fabs(a[r][r]))
                t = i;
        
        for (int j = r; j <= n; j++)
            swap(a[r][j], a[t][j]);
        
        // 归一
        for (int j = n; j >= r; j--)
            a[r][j] /= a[r][r];
        
        // 消为 0
        for (int i = r+1; i < n; i++)
            for (int j = n; j >= r; j--)
                a[i][j] -= a[r][j] * a[i][r];
    }
    
    // 消去上三角
    for (int i = n-1; i >= 0; i--)
        for (int j = i+1; j < n; j++)
            a[i][n] -= a[i][j] * a[j][n];
}

int main() {
    scanf("%d", &n);
    for (int i = 0; i < n+1; i++)
        for (int j = 0; j < n; j++)
            scanf("%lf", &a[i][j]);
    
    // 预处理系数
    for (int i = 0; i < n; i++) {
        double b = 0;
        for (int j = 0; j < n; j++) {
            b += a[i+1][j]*a[i+1][j]-a[i][j]*a[i][j];
            a[i][j] = 2*(a[i+1][j]-a[i][j]);
        }
        a[i][n] = b;
    }
    gauss();
    for (int i = 0; i < n; i++)
        printf("%.3lf ", a[i][n]);
    return 0;
}

开关问题

有 N 个相同的开关，每个开关都与某些开关有着联系，每当你打开或者关闭某个开关的时候，其他的与此开关相关联的开关也会相应地发生变化，即这些相联系的开关的状态如果原来为开就变为关，如果为关就变为开。

你的目标是经过若干次开关操作后使得最后 N 个开关达到一个特定的状态。

对于任意一个开关，最多只能进行一次开关操作。

你的任务是，计算有多少种可以达到指定状态的方法。（不计开关操作的顺序）

题目链接：AcWing 208。

这是一个异或方程组。如果有下列控制：

1
2
3
x1 -> 1, 2, 3
x2 -> 1, 2
x3 -> 1, 3

就有这个方程组，其中 f(x)f(x) 指第 xx 个灯的状态是否改变，方程的解表示每个开关按或不按：

{x1 xor x2 xor x3=f(1)x1 xor x2=f(2)x1 xor x3=f(3)\left\{\begin{aligned} &x_1 \text{ xor } x_2 \text{ xor } x_3=f(1)\\ &x_1 \text{ xor } x_2=f(2)\\ &x_1 \text{ xor } x_3=f(3) \end{aligned}\right.

对于每个自由元，都有两种状态，所以如果方程有解时存在 kk 个自由元，方案是 2k2^k 个。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int N = 35;
int a[N][N], n;

// 返回自由元个数
int gauss() {
    int r, p;
    for (r = 1, p = 1; p <= n; p++) {
        // 找到主元最大的
        int t = r;
        for (int i = r; i <= n; i++)
            if (a[i][p]) {
                t = i;
                break;
            }
        
        if (a[t][p] == 0) continue;
        
        // 交换 a[r], a[t]
        for (int j = 1; j <= n+1; j++)
            swap(a[r][j], a[t][j]);
        
        // 消元
        for (int i = r+1; i <= n; i++)
            if (a[i][p])
                for (int j = n+1; j >= p; j--)
                    a[i][j] ^= a[r][j];
        
        r++;
    }
    
    if (r <= n) {
        for (int i = r; i <= n; i++)
            if (a[i][n+1]) return -1;
        return n-r+1;
    }
    
    return 0;
}

void solve() {
    memset(a, 0, sizeof a);
    scanf("%d", &n);
    
    // a[i][n+1] 表示状态是否发生改变
    // 也就是方程组解的矩阵
    for (int i = 1; i <= n; i++)
        scanf("%d", &a[i][n+1]);
    for (int i = 1; i <= n; i++) {
        int t; scanf("%d", &t);
        a[i][n+1] ^= t;
        // 每个开关一定能控制自己的灯，被这个地方坑了一下
        a[i][i] = 1;
    }
    
    // 第 y 个灯可以被第 x 开关控制
    int x, y;
    while (scanf("%d%d", &x, &y), x || y)
        a[y][x] = 1;
    
    int k = gauss();
    if (k == -1) puts("Oh,it's impossible~!!");
    else printf("%d\n", 1<<k);
}

int main() {
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}