[关闭]
@Kurunie 2015-06-18T19:04:08.000000Z 字数 5672 阅读 4209

浅谈异或方程组——高斯消元

高斯消元 异或方程组

【异或（XOR）运算由于与“奇偶性”密切相关，经常出现在有关奇偶性和二进制的题目中】
很多异或问题都要涉及到解异或方程组，因此首先要搞懂异或方程组的解法。

（1）异或方程组的格式（设有n个未知数，m个方程）：

a11*x1 XOR a12*x2 XOR ... XOR a1n*xn = b1
a21*x1 XOR a22*x2 XOR ... XOR a2n*xn = b2
......
am1*x1 XOR am2*x2 XOR ... XOR amn*xn = bm

其中的所有a、b、x的值均为0或1。解异或方程组就是给出了所有的系数a和b之后，求出解（x1, x2 ... xn）。一般来说，题目的要求无非就是如下几种：

  1. 求任意一组解
  2. 求解的总数
  3. 求出最优解（比如字典序最小的解或者加权以后权值最大/小的解等）

（2）求解异或方程组的离线算法：

整体思想与普通方程组的高斯消元解法类似。
第一阶段（消元）：
设置一个变量p，一开始p=1，然后，i从1到n，每次执行：

  1. 若在第p个及以后的方程中，至少有一个方程的xi系数为1，设为第q个方程，则先将第p、q个方程交换，然后用第p个方程去XOR后面剩下的所有xi系数为1的方程（各系数包括b都对应进行XOR运算，实际上就是矩阵初等行变换），将它们的xi系数均变成0，最后p加1

  2. 否则（第p个及以后的方程中xi系数均为0），xi是自由元（xi不管是取0还是1，都能导出整个方程组的一个解），按照自由元处理（随便定值，然后将前面所有xi系数为1的方程全部代入xi的值XOR掉），p值不变。(感谢HN SYJ神犇，在总结中讲到了如果遇到自由元的话p值是不能变的，否则会疵，见这里）

第二阶段（回代求解）：
在第一阶段结束后，若出现了S个自由元，则会在整个方程组的最后出现S个多余方程，它们的等号左边所有系数均为0（显然值也为0）。对于这些多余方程，如果有等号右边为1的，则整个方程组无解。否则，如果所有多余方程的等号右边均为0，则整个方程组有2S个解（因为S个自由元可以随便取值，都能导出解）。如果要求具体的解，则从第(m-S)个方程开始，往回代，设立变量p，初始p=m-S，逆序枚举每个未知数（自由元除外），对于未知数xi，可以从第p个方程中直接得到其值（因为此时第p个方程等号左边必然只有xi系数为1，等号右边是几，xi值就是几），然后，将之前的所有xi系数为1的方程全部用第p个方程XOR掉，p--。
具体实现的时候有一些技巧：
1. 在第一阶段第<1>步XOR的时候，实际上只要枚举xi及其后面的未知数就行了，因为前面的未知数已经消掉了；在第<2>步XOR的时候，只需要xi和等号右边的系数b，因为其它的系数均为0；
2. 第一阶段的p可以直接留到第二阶段继续用，只不过要减1。

（3）求解异或方程组的在线算法：

这里的在线算法其实指的是，它可以随时在后面增加新的方程，直接扩展已求出的解以及更新解的总数，而不用每次重新求解。这个算法在某些求”XOR值最大“的问题中显然比较有用。
为了讲清楚这个算法，首先要引入”关键元“：在消元过程中，每个方程有至多一个”关键元“，它决定了该方程可以去XOR后面的哪种方程。设C[i]为关键元为xi的方程编号，若C[i]==-1（以xi为关键元的方程不存在），则xi为自由元，反之亦然。一开始，所有的C值均为-1。
在新加入一个方程时，i从1到n枚举，若该方程xi系数为1，则看一下C[i]是否为-1，若为-1，则xi为该方程关键元，C[i]设为该方程编号，结束（此时，xi不再是自由元，因此整个方程组解的总数减半）；若不为-1，则用编号C[i]的方程去XOR该方程（其实只要XOR xi及其之后的部分），将该方程中的xi消去。如果最终i枚举到n时仍然木有找到关键元，则该方程无关键元，也就是说该方程是一个多余方程（0=0或0=1），如果是0=0，则不影响解的总数，如果是0=1，则整个方程组无解了。
如果要求具体的解，则与离线算法的第二阶段类似，只是回代求解时，p不能再是逆序的了，而应该是每个C[i]的值。此外，对于C[i]==-1的（即xi是自由元），随便定值，并代到所有方程中XOR掉它。

显然，以上两种算法的时间复杂度均为O(n2m)。

（4）算法常数优化——压位：

在具体实现的时候，需要用一个bool矩阵A[m][n+1]来存储所有的系数（a、b），这样很浪费，因为完全可以用一个int矩阵，一下存储32位。这种压位思想可以带来常数上的优化，因为在用一个方程去XOR另一个时，时间将变为原来的1/32。
实现起来比不压位的要麻烦一些。具体来说，首先，调用原矩阵的第i行第j列是否为1需要写成A[i][j/32] & (1 << (j % 32))（原矩阵第i行第j列为0当且仅当该值为0），对于等号右边的b系数，在原矩阵中为第n列，因此写成A[i][n/32] & (1 << (n % 32))。在实现过程中，为了避免重复计算，可以设n0=n/32，q=n%32，q0=j/32，q1=j%32，n0和q预先算出，在枚举j（未知数编号）时维护q0、q1。
另外需要严重注意的是：在求出了某个未知数xi的值之后代入到其它方程中时，如果是单个代入不整体XOR，则i和n列都要处理，但如果是整体XOR，则要特判一下i和n列在压位之后是否压到同一位里了，若压到同一位里则只能XOR一次！！！！！！！！！
（其实是可以压64位的，但由于long long的常数较大且在求1<

例题：JSOI2012 arc
建模方法：设一组解(x1, x2 ... xn)为：xi==0当且仅当i在上游。对于原图中给出的有向图，若i的出度为偶数，则需要满足i指向的那些点的x值XOR之和为0即可（不管xi是0还是1）；若i的出度为奇数，则需要满足i及其指向的那些点的x值的XOR之和为1即可（不管xi是0还是1）。这样就转化为求异或方程组的特解问题，假设在求解过程中，对于所有的自由元均定为0。
代码（均为压位之后的）：
离线算法：

  1. #include <iostream>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <time.h>
  6. using namespace std;
  7. #define re(i, n) for (int i=0; i<n; i++)
  8. #define re1(i, n) for (int i=1; i<=n; i++)
  9. #define re2(i, l, r) for (int i=l; i<r; i++)
  10. #define re3(i, l, r) for (int i=l; i<=r; i++)
  11. #define rre(i, n) for (int i=n-1; i>=0; i--)
  12. #define rre1(i, n) for (int i=n; i>0; i--)
  13. #define rre2(i, r, l) for (int i=r-1; i>=l; i--)
  14. #define rre3(i, r, l) for (int i=r; i>=l; i--)
  15. #define ll long long
  16. const int MAXN = 2005, KS = 32, MAXN0 = (MAXN + 1) / KS + 1, INF = ~0U >> 2;
  17. int n, n0, q, A[MAXN][MAXN0];
  18. bool res[MAXN], res_ex = 1;
  19. void init()
  20. {
  21.     freopen("arc.in", "r", stdin);
  22.     int q0 = 0, q1 = 0, x, y;
  23.     scanf("%d", &n); n0 = n / KS; q = n % KS;
  24.     re(i, n) {
  25.         scanf("%d", &x);
  26.         if (x & 1) {A[i][q0] |= 1 << q1; A[i][n0] |= 1 << q;}
  27.         re(j, x) {scanf("%d", &y); y--; A[i][y / KS] |= 1 << (y % KS);}
  28.         if (q1 == KS - 1) {q0++; q1 = 0;} else q1++;
  29.     }
  30.     fclose(stdin);
  31. }
  32. void solve()
  33. {
  34.     int q0 = 0, q1 = 0, p = 0, x;
  35.     ll tmp;
  36.     re(i, n) {
  37.         x = -1;
  38.         re2(j, p, n) if (A[j][q0] & (1 << q1)) {x = j; break;}
  39.         if (x >= 0) {
  40.             if (x != p) {re3(k, q0, n0) {tmp = A[x][k]; A[x][k] = A[p][k]; A[p][k] = tmp;}}
  41.             re2(j, p+1, n) if (A[j][q0] & (1 << q1)) re3(k, q0, n0) A[j][k] ^= A[p][k];
  42.             p++;
  43.         } else {
  44.             res[i] = 1;
  45.             re(j, p) if (A[j][q0] & (1 << q1)) {A[j][q0] &= ~(1 << q1); A[j][n0] ^= 1 << q;}
  46.         }
  47.         if (q1 == KS - 1) {q0++; q1 = 0;} else q1++;
  48.     }
  49.     re2(i, p, n) if (A[i][n0] & (1 << q)) {res_ex = 0; return;} p--;
  50.     rre(i, n) {
  51.         if (q1) q1--; else {q0--; q1 = KS - 1;}
  52.         if (!res[i]) {
  53.             res[i] = A[p][n0] & (1 << q);
  54.             re(j, p) if (A[j][q0] & (1 << q1)) {
  55.                 A[j][q0] ^= A[p][q0]; if (q0 < n0) A[j][n0] ^= A[p][n0];
  56.             }
  57.             p--;
  58.         }
  59.     }
  60. }
  61. void pri()
  62. {
  63.     freopen("arc.out", "w", stdout);
  64.     if (res_ex) {
  65.         int sum = 0; bool SPC = 0;
  66.         re(i, n) if (!res[i]) sum++;
  67.         printf("%d\n", sum);
  68.         re(i, n) if (!res[i]) {
  69.             if (SPC) putchar(' '); else SPC = 1;
  70.             printf("%d", i + 1);
  71.         }
  72.         puts("");
  73.     } else puts("Impossible");
  74.     fclose(stdout);
  75. }
  76. int main()
  77. {
  78.     init();
  79.     solve();
  80.     pri();
  81.     return 0;
  82. }

在线算法：

  1. #include <iostream>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. using namespace std;
  6. #define re(i, n) for (int i=0; i<n; i++)
  7. #define re1(i, n) for (int i=1; i<=n; i++)
  8. #define re2(i, l, r) for (int i=l; i<r; i++)
  9. #define re3(i, l, r) for (int i=l; i<=r; i++)
  10. #define rre(i, n) for (int i=n-1; i>=0; i--)
  11. #define rre1(i, n) for (int i=n; i>0; i--)
  12. #define rre2(i, r, l) for (int i=r-1; i>=l; i--)
  13. #define rre3(i, r, l) for (int i=r; i>=l; i--)
  14. #define ll long long
  15. const int MAXN = 2005, KS = 32, MAXN0 = (MAXN + 1) / KS + 1, INF = ~0U >> 2;
  16. int n, n0, q, A[MAXN][MAXN0], B[MAXN];
  17. bool C[MAXN], res[MAXN], res_ex = 1;
  18. void init()
  19. {
  20.     freopen("arc.in", "r", stdin);
  21.     int q0 = 0, q1 = 0, x, y, _q0, _q1;
  22.     scanf("%d", &n); n0 = n / KS; q = n % KS;
  23.     re(i, n) {
  24.         scanf("%d", &x);
  25.         if (x & 1) {A[i][q0] |= 1 << q1; A[i][n0] |= 1 << q;}
  26.         re(j, x) {
  27.             scanf("%d", &y); y--; _q0 = y / KS; _q1 = y % KS;
  28.             A[i][_q0] |= 1 << _q1;
  29.         }
  30.         if (q1 == KS - 1) {q0++; q1 = 0;} else q1++;
  31.     }
  32.     fclose(stdin);
  33. }
  34. void solve()
  35. {
  36.     re(i, n) B[i] = -1;
  37.     int q0, q1, x;
  38.     re(i, n) {
  39.         q0 = q1 = 0;
  40.         re(j, n) {
  41.             if (A[i][q0] & (1 << q1)) {
  42.                 x = B[j];
  43.                 if (x == -1) {B[j] = i; break;} else re3(k, q0, n0) A[i][k] ^= A[x][k];
  44.             }
  45.             if (q1 == KS - 1) {q0++; q1 = 0;} else q1++;
  46.         }
  47.     }
  48.     q0 = n0; q1 = q;
  49.     rre(i, n) {
  50.         if (q1) q1--; else {q0--; q1 = KS - 1;}
  51.         if ((x = B[i]) >= 0) {
  52.             res[i] = A[x][n0] & (1 << q);
  53.             re(j, n) if (j != x && (A[j][q0] & (1 << q1))) {
  54.                 A[j][q0] ^= A[x][q0]; if (q0 < n0) A[j][n0] ^= A[x][n0];
  55.             }
  56.         } else {
  57.             res[i] = 1;
  58.             re(j, n) if (A[j][q0] & (1 << q1)) {
  59.                 A[j][q0] &= ~(1 << q1); A[j][n0] ^= 1 << q;
  60.             }
  61.         }
  62.     }
  63.     re(i, n) if ((x = B[i]) >= 0) C[x] = 1;
  64.     re(i, n) if (!C[i] && (A[i][n0] & (1 << q))) {res_ex = 0; return;}
  65. }
  66. void pri()
  67. {
  68.     freopen("arc.out", "w", stdout);
  69.     if (res_ex) {
  70.         int sum = 0; bool SPC = 0; re(i, n) if (!res[i]) sum++;
  71.         printf("%d\n", sum);
  72.         re(i, n) if (!res[i]) {
  73.             if (SPC) putchar(' '); else SPC = 1;
  74.             printf("%d", i + 1);
  75.         }
  76.         puts("");
  77.     } else puts("Impossible");
  78.     fclose(stdout);
  79. }
  80. int main()
  81. {
  82.     init();
  83.     solve();
  84.     pri();
  85.     return 0;
  86. }

原文：http://www.cppblog.com/MatoNo1/archive/2012/05/20/175404.html

添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注