解线性方程组——高斯消元の板子

ATP记得它在很久以前看过一点点高斯消元的东西然后做过一点点题目。。但是当时实在是太zz了所以本来就没有很懂这个东西现在更是忘得差不多了。。
所以现在就当重新学一遍了QwQ

一点口胡的解释

高斯消元。。听起来这个名字很高大上实际上它也确实很高大上但是它的操作过程实际上和数学课学的加减消元法是一个道理的啦。

它可以用来求解线性方程组 Ax=b Ax=b，其中 A A是一个n∗nn*n的矩阵， x x和bb是 n∗1 n*1的矩阵。实际上也就是求解由n个未知量和n个方程组成的n元一次方程组，其中A是系数矩阵，A的第i行代表第i个方程，第j列代表第j个未知数对应的系数。而b是常数项，第i行代表第i个方程等号后面的常数项。

高斯消元还可以求解可逆矩阵的逆矩阵。还可以求解矩阵的秩。（秩是啥？我不会= =）

高斯消元的操作过程实际上非常暴力。。总的来说分两个步骤，首先要通过矩阵行变换也就是加减消元消来消去把这些方程变成某个等价形式。这个等价形式的特征体现在系数矩阵A里就是A变成了一个“上三角矩阵”，也就是只有左上-右下对角线及其上方有数字，这条对角线的下面都是0，大概就是这个样的东西：

⎡⎣⎢⎢⎢⎢⎢A1,1000A1,2A2,2⋮0⋯⋯⋱⋯A1,nA2,n⋮An,n⎤⎦⎥⎥⎥⎥⎥

\begin{bmatrix} A_{1,1} & A_{1,2} & \cdots & A_{1,n} \\ 0 & A_{2,2} & \cdots &A_{2,n} \\ 0 & \vdots & \ddots &\vdots \\ 0 & 0 & \cdots & A_{n,n} \end{bmatrix}
这样的话可以看出来第i个方程有n-i+1个系数，那么最后一个方程只有一个系数，也就是经过这样消元以后最后一个方程变成了 An,n∗xn=bn A_{n,n}*x_n=b_n的形式，那么就可以直接得到 xn x_n的值了。得到 xn x_n以后就可以往回代，因为倒数第2个方程有2个系数，其中一个系数是关于 xn x_n的，把得到的 xn x_n的值带回去以后就可以把它消成只有一个未知量的方程，那么就可以得到 xn−1 x_{n-1}。同理把已经求出的未知量都往回代，就可以求出所有的未知数。

贴一个板子。出自BZOJ1013球形空间产生器，是高斯消元的一个板子题：


void Gauss_Eli(){int num;for (int i=1;i<n;i++){num=i;for (int j=i+1;j<=n;j++)if (fabs(A[j][i])>fabs(A[num][i]))num=j;for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);swap(b[num],b[i]);for (int j=i+1;j<=n;j++){if (fabs(A[j][i])<=eps) continue;double t=A[j][i]/A[i][i];for (int k=1;k<=n;k++)A[j][k]-=A[i][k]*t;b[j]-=b[i]*t;}}for (int i=n;i>=1;i--){ans[i]=b[i]/A[i][i];for (int j=i;j>=1;j--)b[j]-=A[j][i]*ans[i];}
}

可以看出高斯消元具体的实现复杂度是 O(n3) O(n^3)的。这里以及下面我们都把将要留在上三角矩阵的第i行的方程称为未知数i的关键方程，因为如果未知数 xi x_i的值是唯一的，上三角矩阵的第i行必须有未知数 xi x_i的系数。

上面的代码大体来说就是先枚举每一个变量i找到它的关键方程，然后把它的关键方程下面所有方程 xi x_i的系数都消成0让它符合上三角矩阵的形式，最后再进行回代。

首先我们要从上到下枚举每一行。设当前枚举到的是第i行，那么前i-1行已经消成了上三角矩阵的一部分，不能再动了。我们要找第i个未知量的关键方程只能从后面找。上面代码的第4..9行干了一个看起来比较奇怪的事情：从后面的方程中找了一个未知数i的系数最大的方程然后换到第i个位置。据说找一个最大的能减少浮点误差？然而比较实质性的作用还是保证选中的这个“准”关键方程里 xi x_i的系数不是0吧。。如果它题目里给的方程排列比较凑巧，第i行上 xi x_i的系数正好是0，那么它肯定不能作为 xi x_i的关键方程，肯定要找一个系数不是零的搞过来用。所以这一步还是比较重要的。

找到了“准”关键方程以后我们要把它变成关键方程。要干的事情就是让它下面的所有方程 xi x_i的系数都变成0。那么就要用到加减消元法了。设第i行下面的某一行为j，如果第j行中 xi x_i的系数不为0，那么我们就要把它变成0。联想一下数学课上解二元一次方程组的方法，我们要把某一个方程中某一个未知量消掉，方法就是把这个方程扩大一个特定的倍数，让要消掉的这个未知量的系数和另一个方程的对应系数相等，再把这两个方程相减就能消掉目标未知数。那么这里目标未知数 xi x_i的系数是 Ai,i A_{i,i}，我们要把它扩大某个倍数让 Ai,i A_{i,i}变成 Aj,i A_{j,i}，然后用第j个方程和第i个方程作差来消掉 xi x_i。显然只需要把第i个方程扩大 Aj,iAi,i A_{j,i} \over {A_{i,i}}倍就能解决问题。那么求出这个系数，将第i个方程包括后面带着的常数项 bi b_i都乘以这个系数然后作差就可以了。

最后一步是回代，回代的时候要从后往前枚举上三角矩阵中的方程，每次求出一个未知数 xi x_i以后都检查前面所有 xi x_i的系数不为0的方程然后把它代进去，相当于进行一次移项，因为 xi x_i求出来以后原本带未知数的项变成了常数项，那么把它移动到方程右边和 b b合并就可以了。最后当这个方程后面的所有方程都处理过了以后，它自己也只带一个未知数了，并且这个未知数的值也可以直接求出来了。

比较简单的题目

1.异或方程组（BZOJ1923、POJ1222、POJ1681、POJ1830、POJ3185、BZOJ1770）

异或方程组是用高斯消元可以解决的一大类问题，实际上就是普通解方程组的变形。在异或运算里面因为当前变量的出现情况只和系数的奇偶性有关，那么可以直接把系数和常数项变成0或1。而加减消元也变成了直接异或，因为不需要计算扩大的倍数，只要异或一下就能把目标未知数消去。

“开关问题”是异或方程组的经典问题，大概就是指定一些灯，每个灯的初始状态或开或关。每个灯有一个开关，每个开关还有一些相关的开关。如果按下一个开关，那么与之相关的开关都会被按下，相连的灯的开关状态都会改变，求让灯达到某个目标状态的方案。

可以发现灯的开关状态只和按下开关次数的奇偶性有关，那么可以构造异或方程组的模型。把每个开关按或没按设成未知量xix_i，如果按了就赋值为1，没按就赋值为0，那么可以发现某个灯状态改没改变取决于它自己的开关和所有与它相连的开关。按照相连状态建立矩阵，对于每个灯建立一个方程，把所有会影响到它的开关对应的系数赋值为1，不会影响到它的就赋值为0。常数项b对应着这个灯的目标状态，然后解异或方程组就可以了。

在解异或方程组的时候还有一个有趣的优化。可以发现在消元的时候实际上就是系数矩阵的每一行按位异或，那么如果用bitset存储系数矩阵A，每次消元的时候就可以一次异或一整行，省去了枚举，可以把时间复杂度从 O(n3) O(n^3)优化到 O(n364) O({n^3 \over {64}})。有些题不用bitset是过不了的。。

这里贴出使用bitset的板子。常数项 bi b_i存在了 Ai,n+1 A_{i,n+1}中。

bitset<250>A[250];
void Gauss_Eli(){int num;for (int i=1;i<=n;i++){if (A[i][i]==0){num=i;for (int j=i+1;j<=n;j++)if (A[j][i]!=0){num=j;break;}swap(A[num],A[i]);}for (int j=i+1;j<=n;j++)if (A[j][i]!=0) A[j]^=A[i];}for (int i=n;i>=1;i--){ans[i]=A[i][n+1];for (int j=i;j>=1;j--)if (A[j][i]!=0)A[j][n+1]=A[j][n+1]^ans[i];}
}

在异或方程组中常常考察的另一个问题就是自由元的问题。在消元的过程中，可能某个未知数 xi x_i根本就找不到“关键方程”，那么 xi x_i就可以随便取值，取什么值都可以构造出一组合法解。这样的 xi x_i叫做自由元。在异或方程组中每个自由元都有0或1两种取值，那么如果有t个自由元，就可以构造出 2t 2^t组合法解。需要注意的一个问题是如果题目要求某种最优方案，一般来说是需要爆搜所有自由元的。比如要求按下开关最少的方案，不一定自由元中1的个数最少，最后构造出的解1的个数就最少。

然而还有一个问题就是当实际操作中出现“自由元”的时候，当前这一行i中变量 xi x_i的系数为0，而后面可能还有不是0的变量。那么随着回代过程的进行，这些变量会被一一消掉然后和常数项合并。最后枚举到第i行的时候，本来应该有一个 xi x_i还没被消掉留在这里的，但现在因为 xi x_i的系数为0所以没有这个变量了，也就相当于方程的左边是0。如果方程右边最后剩下的常数项也是0，那当然非常好， xi x_i就是一个自由元；但如果最后剩下的常数项不是0，说明当前方程组无解，因为出现了形如 0=1 0=1这样不能满足的方程。

2.同余方程组（POJ2065、POJ2947）

同余方程组又是另一个变形。处理也非常简单，把消元过程中的除法换成逆元就可以了。

实际上在面对处理的数字范围比较小而结果又保证是整数的题目的时候，为了避免浮点误差可以采用找一个大质数取模的方法来避免浮点运算。。比如BZOJ4004当时ATP做的时候被卡精度卡得惨不忍睹然后就改了这种方法。。但是实测好像比直接浮点数运算要慢？

高斯消元的重要应用——线性基

关于线性基等线性代数知识可以参考2014年的集训队论文：匡正非《浅谈线性相关》

线性基是什么？

线性基是一组线性无关的向量，所谓线性无关也就是说这些向量无论如何组合，加上什么系数都凑不出0向量。形式化地，一组n维向量 ai=[x1,x2,...,xn](i=1..k) a_i=[x_1,x_2,...,x_n](i=1..k)线性无关当且仅当不存在一组不全为0的系数 [c1,c2,...,ck] [c_1,c_2,...,c_k]使得 c1a1+c2a2+...+ckak=0 c_1a_1+c_2a_2+...+c_ka_k=0。

线性基和我们平常所说的向量基底实际上是一个东西。考虑由一组n维向量定义的线性空间，显然我们可以用n个线性无关的向量表示出这个空间里任意一个向量，当然也能表示出给定的n个向量以及它们的各种组合。这n个线性无关的向量就叫做这个线性空间的基，也就是基底。

举个例子，在最简单的二维平面上，每一个有序实数对 (xi,yi) (x_i,y_i)定义了一个向量。而由平面向量基本定理可得，不共线的任意两个向量都可以作为一组基底来表示这个平面内的所有向量。而两个向量共线的条件正是 xi=kyi x_i=ky_i，即，存在系数使得 xi−kyi=0 x_i-ky_i=0，它们线性相关。这个情况推广到三维甚至更高维空间里也是成立的。

也就是说，如果给定了一些向量，我们只需要求出它的一组线性基，就可以用线性基的组合来表示所有这些向量以及它们可能组合出的向量。

那么这个东西跟高斯消元有什么关系呢？可以发现高斯消元消的那个矩阵A的每一行实际上都可以看成是一个向量。那么我们可以通过控制消元的过程，使得在消元的时候筛去所有能够被其它向量表示的向量，最后选定每一维对应的基底。

先贴一个代码，以求解异或形式下的线性基为例：

    for (int i=1;i<=n;i++){for (int j=63;j>=0;j--)if (a[i].num>>j&1)if (b[j]==0){b[j]=i;break;}else a[i].num^=a[b[j]].num;}

首先我们从1到n枚举每个向量，这里的每个向量用一个long long类型的数字表示。从高位到低位枚举每一位，如果当前向量的这一位是1，那么说明它这一位的数字还没有被前面任何向量消去，就检查这一位有没有记录对应的基底。如果没有记录的话就把当前这个数字作为对应的基底然后退出循环，否则，如果这一位已经找到了对应的基底，说明当前这个数字不应该有这一位，要把它消掉，也就是和这一位找好的基底异或一下把这一位消掉。

这个过程做完以后数组 b[i] b[i]里存储的就是第i位对应的基底。实际上可以发现这个过程和消上三角矩阵的过程有些相似，都是为每一位选定“关键向量”，并且保证关键向量下面的所有向量都没有这一位。那么可以发现，因为每一位都只有一个“关键向量”，那么如果从求出的线性基里面任意选择两个进行异或运算，设为 ai a_i和 aj a_j，不妨设 i>j i \gt j，那么由求解过程可得 ai a_i的第i位是它最高位的1，而 aj a_j的这一位是0。那么 ai a_i和 aj a_j异或就怎么也得不到0。推广到多个向量也是如此。也就是说，通过上面的过程，我们能够求出一组线性无关的向量，也就是给定向量的线性基。

注意到当这个过程结束以后，有些向量会被异或成0。那么这些向量一定可以被原向量组中其它某些向量线性表示，也就是说原向量组中存在一个组合使得它们的异或值为0。

线性基有什么用？

1.贪心

很多题目为每个向量赋了一个权值，要求你选权值最大的线性无关组。求线性无关组当然就是求线性基，但如何保证权值最大呢？这里我们要做的就是把向量按照权值从大到小排序然后按照这个顺序求解线性基，那么选出的就是权值最大的。好像有一个叫做拟阵的东西能证明这个东西？然而拟阵是啥ATP并不会。。。

例题：
1. BZOJ2460 元素
2. BZOJ4004 装备购买
3. BZOJ3105 新Nim游戏

2.求第k大组合

线性基的一个特点就是它有一定的有序性，那么我们可以利用线性基的这种有序性来解决一些统计问题。常见的问题是给出一个向量组，求这些向量组中的向量异或可以得到的第k大向量。因为如果直接在原向量组里进行操作会有很多出现乱七八糟0向量还有重复向量的问题，我们不妨求出线性基，用线性基来代替原向量进行操作。可以发现当求出的线性基只有一个向量的时候，只能表示出1个数字；如果有两个向量，就可以表示出3个数字……以此类推，如果求出的线性基由n个向量组成，就可以表示出 2n 2^n个不同的数字。至于能不能表示出0这个数字需要特判，只有在求解线性基的过程中出现了某个被异或成0的向量才说明原向量中存在异或和为0的组合。

那么如何求第k大呢？可以发现如果把所有求出的线性基从大到小排好序，一共有 n n个向量。而合法的k显然不超过2n2^n，那么我们可以把k二进制分解，从高位到低位对应着从大到小的线性基。如果这一位有1，就选上这一位对应的基底。这样就能构造出第k大组合。

而还有一个问题就是这样求出的“第k大”显然是指重复数字只算一次的第k大，因为线性基异或出来的数字是没有重复的。而如果要计入重复数字怎么办呢？可以发现在求解线性基的时候有一些向量被“筛”掉了也就是最后被异或成了0，这些被异或成0的向量我们可以把它看做是“第0位”的基底。这些基底因为本身值就是0，那么不管异或上几个都不会对最后结果产生影响，也就是如果数字 N N在不计入重复数字的时候是第k个，而一共有tot个“0”基底，那么NN在计入重复数字的时候就是第 N∗2tot N*2^{tot}个。因为每一个不重复数字都可以异或上任意个数的“0”基底得到相同的数字。

例题：
1. HDU3949 XOR
2. BZOJ2115 Xor
3. BZOJ4269 再见Xor
4. BZOJ2844 albus就是要第一个出场
5. BZOJ2322 梦想封印
6. BZOJ4568 幸运数字

【其实还有好多线性基乱七八糟的题目ATP还没有做。。等到做了以后再来填坑= =】

高斯消元求解矩阵的逆

Update-2017.1.12：昨天ATP突然手贱点开了一道叫做BZOJ3168的题目。。然后就各种用高斯消元乱搞结果WA了N遍然后发现自己的愚蠢做法非常麻烦然后上网一查题解发现人家都是求的什么矩阵的逆。。（°A°）

啥叫矩阵的逆？

根据ATP的脑补。。矩阵的逆就和整数的逆元和分数的倒数是一个意义的。。因为在矩阵里面“1”就是单位矩阵 I I，那么如果存在两个矩阵AA和 B B满足A∗B=IA*B=I，就把 B B叫做AA的逆矩阵，记作 A−1 A^{-1}。

这里有一个需要注意的地方，严格的逆矩阵的定义应该是两个矩阵 A A和BB满足 A∗B=I A*B=I并且 B∗A=I B*A=I的，但可以证明这两个条件满足一个就可以推出另外一个，所以矩阵的“逆”这种关系应该说是相互的，即如果有两个矩阵满足 A∗B=I A*B=I，那么即可以说 B B是AA的逆，又可以说 A A是BB的逆。显然互逆矩阵是可交换的。

如何求矩阵的逆？

以下所有操作提到的矩阵都是方阵，即行数等于列数的矩阵。

如果给定一个A矩阵，我们通过某种乱七八糟的变换把它变成了一个单位矩阵，那么这个效果是不是就和乘上了一个 A−1 A^{-1}是一样的呢？我们暂且把这个过程叫做把A矩阵“除以A”的过程。

那么我们在把A“除以A”的过程中如果同时维护一个初始是单位矩阵的C矩阵，A做什么操作那C就跟着做什么操作，是不是也就相当于把C“除以A”了呢？因为C初始是单位矩阵“1”，那么“除以A”以后就得到了“ 1A 1 \over A”矩阵，也就是 A−1 A^{-1}啦。

更正经一点的说，根据矩阵初等变换的原理，如果矩阵 A A能通过初等行变换变成BB，那么存在可逆矩阵 P P使得P×A=BP \times A=B。如果把 A A和BB拼起来形成一个新的矩阵 H=(A,B) H=(A,B)，用 P P左乘HH矩阵得到的结果实际上就是 (PA,PB) (PA,PB)。那么如果我们把 B B设成EE，用 P P左乘HH得到的结果就是 (PA,P) (PA,P)。又因为我们要求 A A的逆矩阵，那么PA=EPA=E，那么我们只要把 A A变成EE，就能把 E E变成A−1A^{-1}啦。

先贴上代码：

void Gauss_Eli(){int num;long long t;for (int i=1;i<=n;i++) Inv[i][i]=1;for (int i=1;i<=n;i++){num=i;for (int j=i;j<=n;j++)if (Abs(A[j][i])>Abs(A[num][i])) num=j;for (int j=1;j<=n;j++){swap(A[i][j],A[num][j]);swap(Inv[i][j],Inv[num][j]);}t=powww(A[i][i],Mod-2);for (int j=1;j<=n;j++){A[i][j]=(A[i][j]*t)%Mod;Inv[i][j]=(Inv[i][j]*t)%Mod;}for (int j=1;j<=n;j++)if (i!=j&&A[j][i]!=0){t=A[j][i];for (int k=1;k<=n;k++){A[j][k]=(A[j][k]-(A[i][k]*t)%Mod)%Mod;Inv[j][k]=(Inv[j][k]-(Inv[i][k]*t)%Mod)%Mod;}}}
}

执行完这个过程以后，Inv数组里面存储的就是A的逆矩阵。基本思路是通过矩阵变换把A矩阵变成单位矩阵，其中肯定会把A的某些行乘以某些系数，那么就把维护的这个单位矩阵也在同样的行上乘以同样的系数就可以了。

具体过程是这样的：

首先仍然要为每一列 i i选定一个“代表行”交换到第ii行来，让 A[i][i] A[i][i]变成1，第i列的其它地方都变成0。那么第一步就是要找到一个系数不为0的行交换过来。而因为在A矩阵里执行了交换操作，那么维护的Inv矩阵也就是初始是单位矩阵的那个矩阵也要执行相同的交换操作。

然后要保证第i列只有第i行有一个1其它地方都是0。首先要把第i行整体除以 A[i][i] A[i][i]，也就是乘以 A[i][i] A[i][i]的逆元。然后把Inv矩阵的第i行也乘以相同的数字。然后扫描所有的n行，如果发现有一行的第i列不是0，我们就要把它变成0。仍然利用“加减消元法”，因为我们刚才已经把 A[i][i] A[i][i]变成1了，所以只要把第i行乘以 A[j][i] A[j][i]，那么第i个位置的系数就从1变成了 A[j][i] A[j][i]。这个时候只要用第j行去减一下就把第i个位置的系数给消掉了。注意这里要乘一个数的时候一行之内所有数都要一块乘，并且要在Inv矩阵里乘以相同的数。

这样搞完了以后A矩阵就变成了单位矩阵，而Inv矩阵里存储的就是A矩阵的逆了。

看到这里可能会出现一个问题：在消元的过程中不是要为每一列选定一个系数不为0的“代表行”吗？如果剩下的所有行在第i列的系数都是0，也就是找不到这样一个合法的“代表行”怎么办？

这个时候就说明，无论做什么操作都不可能把矩阵A变成单位矩阵，也就是A没有逆。在线性代数中，这种矩阵叫做“奇异矩阵”。有逆的矩阵就叫做“非奇异矩阵”。

高斯消元求解模线性方程组

Update-2017.2.24：本来一开始发这篇博客的时候ATP不想把模线性方程组这东西另开一个分类的。。因为一开始觉得不就是把除法换成逆元就可以了么。。事实证明ATP还是too naive。。。要是没有逆元怎么办？ATP表示两眼抓瞎。。。

一开始ATP在取模意义下是这么写高斯消元的：

void Gauss_Eli(int n){for (int i=1;i<=n;i++){int num=i;for (int j=i+1;j<=n;j++)if (abs(A[j][i])>abs(A[num][i]))num=j;for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);for (int j=i+1;j<=n;j++)if (A[j][i]!=0){long long t=A[j][i]*getinv(A[i][i])%Mod;for (int k=1;k<=n;k++)A[j][k]=(A[j][k]-A[i][k]*t)%Mod;}}
}

显然如果不良心的出题人把Mod换成一个随便的奇怪数字没法求逆元就只能两眼抓瞎了。。如何不抓瞎？关键就是要避免除法。既要避免除法，又要用第i行消掉第j行的 A[i][i] A[i][i]。。

但是我们如果再考虑一下这个东西，发现它能够成功消元的根本原因就是找到了一个“相等量”，在这里就是利用 A[j][i] A[j][i]这个“相等量”。如何利用呢？关键就在于上面计算的系数 t t。因为t=A[j][i]A[i][i]t={{A[j][i]}\over {A[i][i]}}，把这个系数t乘到第 i i行的时候就能成功达到了把A[i][i]A[i][i]转化为 A[j][i] A[j][i]的目的，再用第j行去减一下就可以利用第 i i行消掉第jj行了。

因为不能做除法的限制，我们不能再利用 A[j][i] A[j][i]这个“相等量”了，那么我们就考虑寻找另外一个“相等量”。实际上理论上可行的“相等量”有很多，比如我们只需要无脑地把第i行乘上 A[j][i] A[j][i]，第j行乘上 A[i][i] A[i][i]，就构造出了 A[j][i]×A[i][i] A[j][i] \times A[i][i]这个相等量。

不过好像大家都用的是 lcm(A[j][i],A[i][i]) lcm(A[j][i],A[i][i])这个东西当做“相等量”。。。ATP也不知道这有什么科学依据但是也没试过别的方法。。把第i行和第j行分别乘上 lcmA[i][i] lcm \over A[i][i]和 lcmA[j][i] lcm \over A[j][i]就可以了。。

贴代码。。。。

void Gauss_Eli(int n){int num;for (int i=1;i<=n;i++){num=i;for (int j=i+1;j<=n;j++)if (Abs(A[j][i])>Abs(A[num][i])) num=j;for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);for (int j=i+1;j<=n;j++)if (A[j][i]!=0){long long lcm=get_LCM(A[i][i],A[j][i]),ta,tb;ta=lcm/Abs(A[j][i]);tb=lcm/Abs(A[i][i]);if (A[j][i]*A[i][i]<0) tb=-tb;for (int k=1;k<=n;k++)A[j][k]=(ta*A[j][k]-tb*A[i][k])%Mod;}}
}

这样就完美的避免了除法！妈妈再也不用担心我求不出逆元辣！

高斯消元求解行列式

Update-2017.2.24：某天，ATP手贱看了一个叫做基尔霍夫矩阵的东西。。然后想这玩意儿不就拿高斯消元一消的事儿吗。。然后ATP刷了两道水题感到自我感觉良好。。

第二天ATP点开了一个叫做BZOJ4031的题目。。第一眼：woc这不就是个裸题。第二眼：woc这什么鬼模数出题人又搞事情。。

当时ATP知道的取模意义下做高斯消元的方法只有求逆元这一种。。接下来ATP上网get到了用最小公倍数乱搞的方法，但是发现求出来行列式不对了。。思考N久无果以后继续上网乱扒，终于扒出了另一种十分科学的求行列式的方法。。

那么这里就从最基本的开始，一步一步解释ATP是怎么掉坑里的吧。。

最基本的求解方法——硬上double

高斯消元求解行列式的方法就是把它消成上三角矩阵，然后这个矩阵主对角线的乘积再乘上 (−1)s (-1)^s再取个绝对值—— s s是高斯消元过程中交换两行的次数——就是矩阵的行列式。原因在各种线性代数教材和网上都有严格的证明，ATP这里只是顺便提一下。。

行列式的计算公式是

∑p为1−n的排列(−1)σ(p)∏i=1nai,pi

\sum_{p为1-n的排列} (-1)^{\sigma(p)}\prod_{i=1}^{n}a_{i,p_i}
其中 σ(p) \sigma(p)是排列p的逆序对数目。

可以看出计算过程中有一步就是选出不同行并且不同列的数字乘起来。那么如果我们要求行列式的这个东西是个上三角矩阵，那么如果我们在某一行选择了一个在主对角线上方的元素，一定会在另外一行选一个主对角线下方的元素。而主对角线下方的元素都是0，这一项对答案没有贡献。也就是说只有选择的n个元素都在主对角线上才能对答案产生贡献。

根据行列式的性质，把行列式的任意两行或者两列交换一次，行列式的值会取反。而把行列式的某一行或者某一列乘上一个数加到另外一行去，行列式的值不变。那么用高斯消元就可以直接搞了。对于数据规模比较小的问题直接硬上double就可以。double的另一个便利之处就是根本不需要维护交换次数，因为反正最后要的是绝对值，符号是啥就不用管了。

又要取模？又没有逆元？

就是有些不良心的出题人给你像1e9这样的模数。。加个7多愉快。。

一般要取模的题如果硬上double不取模就直接炸到不知道哪里去了。。但是取模的话就要求逆元。。

什么？我们可以不做除法？但是问题就在于这里要的是行列式的值。。在那个不做除法的方法里我们有一个最根本的缺陷就在于消掉第 j j行的A[j][i]A[j][i]的同时还把 A[j][i] A[j][i]本身给乘上了一个系数。根据行列式的性质，如果行列式的某一行全体乘上了一个实数 t t，那么行列式就会由DD变成 t∗D t*D。要还原出原来的 D D就又跟除法扯上关系了。。。

在这种情况下要解决这个问题，就又出现了一种高斯消元的新姿势：辗转相除法！

先贴代码再解释：

void Gauss_Eli(int n){int num;for (int i=1;i<=n;i++){num=i;for (int j=i+1;j<=n;j++)if (Abs(A[j][i])>Abs(A[num][i])) num=j;if (num!=i) mak^=1;for (int j=1;j<=n;j++) swap(A[i][j],A[num][j]);for (int j=i+1;j<=n;j++)while (A[j][i]!=0){long long tmp=A[j][i]/A[i][i];for (int k=1;k<=n;k++)//做完一次操作以后A[j][i]相当于对A[i][i]取模了A[j][k]=(A[j][k]+Mod-tmp*A[i][k]%Mod)%Mod;if (A[j][i]==0) break;mak^=1;//注意要维护取反标记，不能直接取绝对值for (int k=1;k<=n;k++) swap(A[j][k],A[i][k]);}}
}

仍然回到我们的根本目的：利用A[i][i]A[i][i]消掉 A[j][i] A[j][i]。现在构造相等量的方法已经行不通了，那么我们再退一步：考虑还有什么可行方法能够把 A[j][i] A[j][i]变成0？

这里利用的是欧几里得算法的原理。在用辗转相除求解gcd的时候我们可以发现，对于一对数字 (a,b) (a,b)，如果我们不断地把 a a变成bb，把 b b变成a%ba \% b，那么总能把其中一个数字变成0，现在我们要利用的就是这个性质。我们要维护的一对数字是 (A[j][i],A[i][i],) (A[j][i],A[i][i],)。我们要做的就是不断把 A[j][i] A[j][i]对 A[i][i] A[i][i]取模然后交换它们的位置，这样就把这对数字换成了 (A[i][i],A[j][i]%A[i][i]) (A[i][i],A[j][i]\%A[i][i])。

具体在矩阵里面做的时候我们肯定要同时对第i行和第j行进行操作，操作的根本目的还是对这两行的第i个元素进行变换，那么我们每次要把 A[j][i] A[j][i]变成 A[j][i]%A[i][i] A[j][i]\%A[i][i]，设 ⌊A[j][i]A[i][i]⌋=t \lfloor{A[j][i] \over A[i][i]}\rfloor=t，那么我们要做的就是把第i行乘上t然后用第j行去减。然后交换这两行的位置，一直到 A[j][i] A[j][i]变成0为止。这样的话利用的只是初等变换 rj−t×ri r_j-t\times r_i，还没有用到除法，就能够解决这个问题了。

注意在取模意义下用高斯消元解行列式的时候一定要维护交换次数来得到符号，因为在取模意义下全都是正数，取绝对值是没有意义的，所以要用专门记录的符号位来区分正数和负数。