昨天有同学在别的群里问了一道题：4×4的格子，每个格子 有黑和白两种状态，每次点击一个格子，它和它相邻的所有格子 都会同时发生翻转，问怎么把所有格子恢复到全白。 以前大学时做这道题都是用枚举第一行策略、推导剩余行、验证最后一行的方法做的，虽然可以剪枝减少一些计算量但是其实意义不太大。后来研究了一段时间矩阵（尤其是异或矩阵），一直也没再遇到过这道题，昨天既然被人问了，那就刷一刷吧，翻出了POJ的三道题，大差不差，略有些细节区别，难度从低到高：
- 1222：限定4乘4，可以到全黑或全白，只求一个最优解的操作的数量。
- 1753：限定5乘6，只能到全0，没有要求最优解（但也没有Special Judge，会不会有多解呢？下文会解释这个问题）。
- 3279：W和H小于等于15，只能到全0，求一个字典序最小的最优解。
最优解定义为 操作的次数最少（解中1的个数最少）
这些题用枚举第一行的方法都可以AC，但既然来写文章，我们肯定不会这么水的过去，今天我们来讲讲数学方法解这道题。
翻转问题的基础
在1222的描述中，已经说过了这个问题的一些基本策略，我们先简单回顾一下：
- 操作的顺序和结果无关（最终只看每个格子被改变次数的“奇偶性”）
- 任何一个格子都没有必要操作两次（同上）
在1753和3279的输出中，我们也可以看出，我们要求的也是一个矩阵，矩阵的每一个元代表了这个格子是否需要操作一次。
线性同余方程
考虑我们要求解的目标：每个格子是否需要操作一次，以及它对应的结果：操作完成之后每个格子的状态是否发生了变化。这两个矩阵的关联呢？则是“操作每个格子 影响到了哪些格子”。这个过程其实可以用一个线性同余方程组来表示（我们这里以2×2来举例）：
x0 * a00 + x1 * a01 + x2 * a02 + x3 * a03 = r0 (mod 2)
x0 * a10 + x1 * a11 + x2 * a12 + x3 * a13 = r1 (mod 2)
x0 * a20 + x1 * a21 + x2 * a22 + x3 * a23 = r2 (mod 2)
x0 * a30 + x1 * a31 + x2 * a32 + x3 * a33 = r3 (mod 2)
这里，x0至x3 分别代表2×2的四个格子 是否需要操作一次，r0至r3 分别代表2×2的四个格子是否需要被改变一次状态，aij表示按下第j个格子，是否会改变第i个格子的状态，因此aij矩阵是已知的：
2乘2的矩阵
1 1 1 0
1 1 0 1
1 0 1 1
0 1 1 1
3乘3的矩阵
1 1 0 1 0 0 0 0 0
1 1 1 0 1 0 0 0 0
0 1 1 0 0 1 0 0 0
1 0 0 1 1 0 1 0 0
0 1 0 1 1 1 0 1 0
0 0 1 0 1 1 0 0 1
0 0 0 1 0 0 1 1 0
0 0 0 0 1 0 1 1 1
0 0 0 0 0 1 0 1 1
上面的方程组，我们可以用向量和矩阵来表示(如果你缺乏相应知识，可以去补习线性代数！)：
A * X = R (mod 2) // 如果X和R是列向量，矩阵元的形式和前面的表示一致。
X * A = R (mod 2) // 如果X和R是行向量，和上面这条等价，但是A矩阵需要转置。
这种表示法称作同余矩阵（向量也是一种特殊的矩阵）。特别的，模2的时候，也可以被称作异或矩阵，因为 a xor b = (a+b) mod 2。
同余矩阵的乘法也满足结合律（注意矩阵乘法不满足交换律），同余矩阵在许多方面也和普通矩阵有着相似的性质，因此我们可以利用这一点来做一些奇怪的事情。
矩阵的秩，解的存在性判定
在线性代数中，一个矩阵A的列秩是A的线性独立的纵列的极大数目。类似地，行秩是A的线性无关的横行的极大数目。即如果把矩阵看成一个个行向量或者列向量，秩就是这些行向量或者列向量的秩，也就是极大无关组中所含向量的个数。
矩阵的秩是其列秩和行秩的较小值。方阵的行秩和列秩总是相等。
线性方程组 解的存在性的关联：
- 秩(系数矩阵)<秩(增广矩阵)：不存在解
- 秩(系数矩阵)=秩(增广矩阵)：存在解：
- 秩(系数矩阵)=未知数的个数：总是有唯一解
- 如果未知数的个数=行数，称作“满秩”，此时系数矩阵存在逆矩阵
- 秩(系数矩阵)<未知数的个数：
- 所有的解可以表达为一个特解 和 K个线性独立的自由基 的线性组合，其中K=未知数的个数-秩(系数矩阵)
- 秩(系数矩阵)=未知数的个数：总是有唯一解
增广矩阵表示把方程每一行的结果 作为矩阵的最后一列出现。
这里特别说一下：对于同余矩阵，上述结论未必成立，满秩也未必代表存在逆矩阵，因为线性同余计算过程中乘法运算未必可逆，无法直接定义一个有效的除法运算。譬如 模4的 1×1矩阵 2 就不存在逆矩阵，因为0到3之间没有任何整数乘以2模4等于1。但是幸运的是这个问题对于异或矩阵来说不存在，满秩的异或矩阵总是存在逆矩阵的，异或矩阵解的存在性和秩也符合上述关联。
上面的描述比较抽象，我们可以具体的来看一看（这里用普通线性方程举例）
方程1
x0 * 1 + x1 * 2 = 4
x0 * 4 + x1 * 4 = 8
方程2:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 6
方程3:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 8
秩是什么意思呢，行秩(系数矩阵)代表了有多少个行，其系数不能通过前面的行的线性组合来得到。观察上面的三个方程：
- 方程1的第二行的系数不能通过第一行的线性组合来得到，系数矩阵的秩为2。有两个未知量，方程的解唯一。
- 方程2的第二行可以通过第一行✖️2来得到，因此系数矩阵的秩只有1。但是考虑增广矩阵就无法由第一行得到了，因此增广矩阵的秩仍为2，该线性方程组无解。
- 方程3的第二行不论是系数矩阵还是增广矩阵，都可以由第一行得到，因此系数矩阵和增广矩阵的秩都为1，且小于未知数数量，该线性方程组有无穷多组解。
这是为什么呢？因为：
方程1
x0 * 1 + x1 * 2 = 4
x0 * 4 + x1 * 4 = 8
# 系数不能由前面方程得到，意味着这是一条独立的方程，可以参与求解。
方程2:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 6
# 系数可以由前面方程得到，但增广不行
# 意味着从前面的方程得到x0*2+x1*4=8，6=8矛盾，因此不存在解。
方程3:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 8
# 系数和增广都可以由前面方程得到
# 意味着从前面的方程就可以得到x0*2+x1*4=8，这行方程去掉也不影响解的数量，等价于只有一行方程有意义。
上面的例子中，行秩就等于秩。如果行数大于等于未知数的数量呢？行秩还等于秩吗？前面的结论还成立吗？
高斯消元
上文讲了秩的定义和解的判定，但没有讲怎么求矩阵的秩，也没有讲到底怎么求出解来。事实上，这两件事都可以用同一个算法来求解，称作高斯消元。
先回忆我们求解线性方程组的消元法（还记得初中的N元一次方程吗……，没错就是那个）
方程1
x0 * 1 + x1 * 2 = 4
x0 * 4 + x1 * 4 = 8
第一行乘以2，得到
x0 * 2 + x1 * 4 = 8
第二行减去第三行，得到
x0 * 2 + x1 * 0 = 0
所以求解出x0 = 0，并可继续求解x1
高斯消元本质上就是把这个过程变成一个规范化的过程，它这么规定，首先通过消元把系数矩阵变成一个三角矩阵（即左下方三角区域全0）
方程1
x0 * 1 + x1 * 2 = 4
x0 * 4 + x1 * 4 = 8 <-- 减去第一行*4
x0 * 1 + x1 * 2 = 4
+ x1 * -4 = -8
第二步，把系数变成一个对角矩阵
x0 * 1 + x1 * 2 = 4 < --减去第二行*-0.5
+ x1 * -4 = -8
x0 * 1 + = 0
+ x1 * -4 = -8
第三步，约去系数
x0 = 0
x1 = 2
从而完成方程的求解。
注意在第一步的时候，我们可能会遇到一些奇怪的情况，譬如：
x1 * 2 = 4 <-- 这一行是x0系数为0，怎么办呢？
x0 * 4 + x1 * 4 = 8
只要x0在后续行的系数还有非0值，那就有很多办法可以解决这个问题，最常见的方法有两种：
x0 * 4 + x1 * 6 = 12 <-- 可以加上该行
x0 * 4 + x1 * 4 = 8
x0 * 4 + x1 * 4 = 8 <-- 可以交换这两行
x1 * 2 = 4
如果一个系数列在后续的所有行都为0呢？很不幸，这意味着这个未知量是一个“自由变量”，也就是说：
- 系数矩阵的秩可能会小于未知量的个数，这个自由变量的取值不影响解的存在性
- 如果后续继续进行第二步，就可以将这一行的所有系数都变成0
- 此时常数列的对应行如果非0，就代表方程无解。
- 如果所有“系数全为0的行”常数列也为0，就代表方程可能有无穷多组解
- 如果存在解：
- 一定存在一个“全部自由变量都为0”的特解
- 这个“自由变量”以及它在之前方程中的系数*-1可以组成一个“自由基”，任何一个方程的解加上这个自由基的倍数都一定还是方程的解。
- 特别的，如果每个自由基仅包含一个自由变量，那么这些自由基之间必然是线性无关的，自由基的数量 等于 自由变量的数量 等于 未知量的个数减去矩阵的秩
注意存在自由变量时，消元无法消到对角矩阵，自由变量所在列的上方可能还存在非0系数。这些非0系数与自由基相关。
我们来具体看一下方程2和3的高斯消元结果：
方程2:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 6 <-- 减去第一行*2
x0 * 1 + x1 * 2 = 4
x0 * 0 + x1 * 0 = -2 <-- 所有系数为0，第二步不涉及此列。
# 常数列非0，代表方程无解。
方程3:
x0 * 1 + x1 * 2 = 4
x0 * 2 + x1 * 4 = 8 <-- 减去第一行*2
x0 * 1 + x1 * 2 = 4
x0 * 0 + x1 * 0 = 0 <-- 所有系数为0，第二步不涉及此列。
# 常数列为0，方程有无穷多组解，x1可以作为自由变量。
# 一定存在一组x1=0的特解，此时x0=4,x1=0
# 该列系数*-1 加上自由变量所在列为1 一组自由基，即为[-2, 1]，意味着如果x0,x1是方程的解，对于任意的c，x0+(-2)*c, x1+1*c都是方程的解（称作通解）
# 该方程只有这一组自由基，意味着它的任何一个解都必须符合4+(-2)*c, 0+1*c
注意有的高斯消元算法的描述会消去这一行，得到的矩阵并非对角矩阵。这更利于笔算但不利于计算机计算，因此我们这里采用保留这一行，且对应行系数全部为0的做法。这两种做法是等价的。
对于满阶方阵来说，高斯消元法还可以用于求解方阵的逆矩阵。我们首先在系数矩阵右侧补充一个单位矩阵；然后可以注意到高斯消元的任何一部操作（不论是交换、线性叠加、约去系数）都等价于左乘一个矩阵 Pi，并且满秩矩阵高斯消元后可以得到一个单位矩阵（对角线系数为1，其余矩阵元全为0）因此操作过程相当于：
Pn*(...(P2*(P1* ([A, I]))))
= Pn * ... * P1 * ([A,I])
= [I, B]
设Pn*...*P1 = V，因此有：
V * [A, I] = [I, B]
可以得到：
V * A = I
V * I = B => V = B
因此B * A = I，B和A互为逆矩阵。
即使矩阵不满秩，也就是不存在逆矩阵，我们在消元过程中得到的“消元矩阵”仍然可以用于快速对一个新的问题消元，从而可以通过预计算降低每次判定的时间复杂度。可惜的是POJ的这几道题没有这样的应用机会。
上述过程和结论对于异或矩阵也都基本符合（但对于非质数的线性同余方阵，因为没有合适的除法运算，所以可能无法进行消元，也就无法应用高斯消元法了）。
但是对于异或矩阵来说，自由基的系数也只能是0或1，因此不是无穷多组解，而是pow(2, 自由变量数量)组解。
解题
在编程表示异或矩阵的时候，我们通常可以用一整个int来表示矩阵的一行（如果int位数不够的话，用多个int，或者使用STL的std::bitset 也可以！），这样我们可以用 AND运算 来快速运算 按位的乘法，用 XOR运算 来快速运算按位的异或。这使得异或矩阵的编程可以极其简练。
并且异或矩阵的取值只有0和1，因此不必进行高斯消元的第三步（约去系数）
// 注意这里假设unsigned int足够表示一行。如果不足够的话，需要将行封装到多个unsigned int里并依次进行异或
// 高斯消元的第一步
for (int i = 0; i < N; i++) {
unsigned int mask = (1U<<i);
if ((A[i] & mask) == 0) {
for (int j = i+1; j<N; j++) {
if (A[j] & mask) {
// 这里采用线性叠加上一个非0行的方法，因为交换似乎开销更大
A[i] ^= A[j];
B[i] ^= B[j];
break;
}
}
}
if ((A[i] & mask) != 0) {
for (int j =i+1; j < N;j++) {
if (A[j] & mask) {
// 消去后面行的该列
A[j] ^= A[i];
B[j] ^= B[i];
}
}
}
}
// 第二步
for (int i = N-1; i>0; i--) {
unsigned int mask = (1U<<i);
for (int j = 0; j < i; j++) {
if (A[j] & mask) {
// 仅当非自由时，消去前面行的该列。
A[j] ^= A[i];
B[j] ^= B[i];
}
}
}
消元之后，我们可能会得到两种结果，或者得到一个单位矩阵（对角全为1），那么满秩&存在逆矩阵，或者部分行系数为0，那么不满秩（可能无解或有多组解）
具体到题目：
- 1753: 有4个自由基，意味着要么无解 要么有16种解，自由基全为0的特解可以根据输入得到（可能无解），然后根据自由基遍历出所有的解并选择最佳的一个。可以切换到白或切换到黑，也可以理解为把问题本身白黑对调再算一遍(和0xffff异或)，取最优解即可。
- 1222：矩阵恰好满秩，因此可以直接用高斯消元法计算出逆矩阵，直接将输入乘以逆矩阵得到解。逆矩阵与问题输入无关，可以预先计算好。
- 3279：题目的大小不确定，意味着满秩和不满秩都可能会遇到，也无法预先计算。因此统一采用高斯消元法直接求解，得到所有的自由基遍历出所有的解并选择最佳的一个。
