用定义求矩阵的行列式
在学习《线性代数》课程的过程中，在面对求给定方阵的行列式的问题时，我们往往采用「高斯消元法」将行列式化为「上三角形行列式」，然后将对角线上的诸元相乘来得到其值。这种方法便于机器计算，但也会不可避免地遇到一个问题——浮点数精度问题。伴随着浮点类型的使用，「高斯消元法」每消去一行，都会对其他有效数字造成误差。因此，本文做出一次「返璞归真」的尝试，在 C 语言中用定义计算方阵的行列式。
所谓「定义」
所谓「定义」，指的是《线性代数》教科书上对「\(n\) 阶行列式」的概念进行叙述时所引入的公式：
\[\det\left[\begin{matrix}a_{11}&a_{12}&\cdots&a_{1n} \\ a_{21}&a_{22}&\cdots&a_{2n} \\\vdots&\vdots&&\vdots \\ a_{n1}&a_{n2}&\cdots&a_{nn}\end{matrix}\right]=\sum(-1)^{t(p_1p_2\cdots p_n)}a_{1p_1}a_{2p_2}\cdots a_{np_n}\]
其中，\(p_1p_2\cdots p_n\) 是 \(1, 2, \cdots, n\) 的排列，\(t(p_1p_2\cdots p_n)\) 是这个排列的逆序数。
使用这条公式来计算行列式，对于使用纸笔进行的计算，这不是一个好方法。但对于我们的机器来说，这种方法或许是值得考虑的。为了让机器能以这种方法解决这个问题，我们将这个问题分解成两大部分：
首先，我们需要生成 \(1, 2, \cdots, n\) 的全排列。仔细分析上面的公式，我们不难发现，等式的右侧正是对 \(1, 2, \cdots, n\) 的每一种排列进行计算和求和。因此，我们需要通过某种手段，生成指定 \(n\) 的自然数列的全排列。在本文中，我们利用递归的方法实现这个部分。
接着，我们需要计算给定数列的逆序数。这对应着 \(t\) 函数——求排列逆序数的函数。这部分我们可以选择「死算」，也可以选择「巧算」——在本文中，我们使用归并排序的思路实现这个部分。
然后，我们进行计算求和。
下面，我们逐一实现这几个部分。
全排列的产生
这部分我们可以抽取出问题：给定一个数组 arr[n]，希望在代码的合理位置得到它的每一个排列。对于这个问题，我们可以使用如下的思路：求 \(n\) 个数的全排列，可以分解成求「第一个数与后面每一个数交换之后，后面 \((n-1)\) 个数的全排列」。具体而言，过程如下：
固定住
arr[0]，然后逐个遍历arr[0]到arr[n - 1]的每一项。假设遍历到第 \(i\) 项，将这一项与
arr[0]交换数值。对于从
arr[1]开始，到arr[n - 1]结束的，长度为 \((n-1)\) 的数组，重复使用本过程，计算这个子列的全排列。即：固定住arr[1]，从arr[1]遍历到arr[n - 1]……直到固定的数字到达了数组的末端。这时，得到一种排列。再次交换第 \(i\) 项和
arr[0]。
以数列 \(1, 2, 3\) 为例，我们演示一下这个过程以帮助理解。（请调整页面显示的宽度，以便下方的过程示意图能完整显示）
“b”代表“固定住”的值。
“^”代表遍历到的值。
水平方向上不同列表示不同层次的调用。
(开始)
1 2 3
b^
(交换)
1 2 3 ===> 1 2 3
b^ b^
(交换)
1 2 3 ===> 1 2 3 (已经遍历到尾，输出 1)
b^ b^
1 2 3 <========
b^
(换回)
1 2 3
b^
(遍历下一项)
1 2 3
b ^
(交换)
1 3 2 ===> 1 3 2 (已经遍历到尾，输出 2)
b ^ b^
1 3 2 <========
b ^
(换回)
1 2 3
b ^
1 2 3 <==== (已经遍历到尾，结束此轮)
b^
(换回)
1 2 3
b^
(下一项)
1 2 3
b ^
(交换)
2 1 3 ===> 2 1 3
b ^ b^
(交换)
2 1 3 ===> 2 1 3 (已经遍历到尾，输出 3)
b^ b^
2 1 3 <=========
b^
(换回)
2 1 3
b^
(下一项)
2 1 3
b ^
(交换)
2 3 1 ===> 2 3 1 (已经遍历到尾，输出 4)
b ^ b^
2 3 1 <========
b ^
(换回)
2 1 3
b ^
2 1 3 <==== (已经遍历到尾，结束此轮)
b ^
(换回)
1 2 3
b ^
(下一项)
1 2 3
b ^
(交换)
3 2 1 ===> 3 2 1
b ^ b^
(交换)
3 2 1 ===> 3 2 1 (已经遍历到尾，输出 5)
b^ b^
3 2 1 <========
b^
(换回)
3 2 1
b^
(下一项)
3 2 1
b ^
(交换)
3 1 2 ===> 3 1 2 (已经遍历到尾，输出 6)
b ^ b^
3 1 2 <========
b ^
(换回)
3 2 1
b ^
3 2 1 <==== (已经遍历到尾，结束此轮)
b ^
(换回)
1 2 3
b ^
(已经遍历到尾，结束此轮)
(结束)
显然，面对这种存在「调用自身过程」的问题，我们的第一反应就是使用递归。我们可以用指针 pHead 和 pTail 来分别指向调用的起点（上面过程中 b 指向的数）和终点（上面过程中数列的尾部）。每次重复调用时，在 pHead 的位置传入 pHead + 1；而终止函数的条件，则是 pHead 与 pTail 相等。按照上面的分析，不难写出如下的函数：
void makePermutation(int* numList, int pHead, int pTail)
{
if (pHead == pTail)
{
// 此时的 numList 数组便是一种排列了。
// 做我们需要做的事。
return;
}
int tempExchange;
for (int i = pHead; i <= pTail; i++)
{
tempExchange = numList[pHead];
numList[pHead] = numList[i];
numList[i] = tempExchange;
makePermutation(numList, pHead + 1, pTail);
tempExchange = numList[pHead];
numList[pHead] = numList[i];
numList[i] = tempExchange;
}
}
可以对上述函数进行完善后验证算法正确性。
求逆序数
哈工大《线性代数》课本上对「逆序数」的定义是：对于排列 \(p_1p_2\cdots p_n\) ，我们把排在 \(p_i\) 前面且比 \(p_i\) 大的数的个数 \(t_i\) 称为全排列 \(p_1p_2\cdots p_n\) 中 \(p_i\) 的逆序数，把这个排列中各数的逆序数之和称为这个排列的逆序数，记为 \(t(p_1p_2\cdots p_n)\)。
看到这个定义的一瞬间，我们首先想到的是利用两个 for 循环嵌套来求取逆序数。这样显然是一种可行的方案；然而，两个 for 循环嵌套这种结构，往往意味着更多的时间开销——毕竟，时间复杂度为 \(O(n^2)\) 呢。
那么有没有别的方法来求逆序数呢？答案是有的。我们使用「归并排序」的思路来求逆序数。
「归并排序」是一种较为高效且具有稳定性的排序方式，它的时间复杂度为 \(O(n\log n)\)。它的思路仍然是一种「二分」法：将数列从中间分开，分成左右两边；两边进行排序之后，将两边的数列「归并」成一个有序的大数列。而对于两边数列的排序，则依然采用归并排序法——显然，这也是一种递归的思路。这篇博客以图解的形式形象地表现了这个过程。
我们明明要求的是逆序数，这和「归并排序」有什么关系呢？我们不妨看下面这个例子：
我们考虑数组 1 2 6 7 3 5 7 8，这个数组的“前半段”为 1 2 3 6，后半段为 3 5 7 8，这两个半段都已经排好序了。现在我们将其进行归并。
1 2 6 7 3 5 7 8
^i ^j
因为 1 < 3，所以选用 i 对应的数：1
1 2 6 7 3 5 7 8
^i ^j
因为 2 < 3，所以选用 i 对应的数：1 2
1 2 6 7 3 5 7 8
^i ^j
因为 6 > 3，所以选用 j 对应的数：1 2 3
1 2 6 7 3 5 7 8
^i ^j
因为 6 > 5，所以选用 j 对应的数：1 2 3 5
1 2 6 7 3 5 7 8
^i ^j
因为 6 < 7，所以选用 i 对应的数：1 2 3 5 6
1 2 6 7 3 5 7 8
^i ^j
因为 7 = 7，所以选用 i 对应的数：1 2 3 5 6 7
i 已经指向前半段的尾部，将后半段从 j 开始的元素直接复制，得到答案：
1 2 3 5 6 7 7 8
我们重点关注这个过程中，「选用 j 对应的数」的步骤。之所以「选用 j 对应的数」，是因为 i 指向的数比 j 指向的数要大。前面的数比后面的数要大……这不正与「逆序数」有某种关联吗？正是如此。
由于我们已经默认归并操作之前，前后两段数列都是顺序的，因此 i 指向的每个数的逆序数都是 0，换言之，逆序数只会由 j 指向的数来「贡献」。而如果 i 指向的数小于等于 j 指向的数，这对于 j 指向的数来说也不能贡献逆序数。只有当 i 指向的数大于 j 指向的数时，对于 j 指向的数才有逆序数的贡献，而这个贡献是多少呢？由于前半段数列已经是顺序的了，既然 i 指向的数已经大于了 j 指向的数，那么前半段数列中 i 之后的所有数都是大于 j 指向的数的。因此这个贡献的值，应该是从 i 开始到前半段数列的尾这一小段数列中数字的个数。
看晕了么？还是以上面的例子来看。
1 2 6 7 3 5 7 8
^i ^j
因为 6 > 3，所以选用 j 对应的数：1 2 3
这一步，我们选用了 j 对应的数：3。那么对于这个 3 而言，它的逆序数则由前半段数列中，从 i 指向的数开始的所有数—— 6 和 7 ——贡献，为 2。
1 2 6 7 3 5 7 8
^i ^j
因为 6 > 5，所以选用 j 对应的数：1 2 3 5
这一步，我们选用了 j 对应的数 5。那么对于这个 5 而言，它的逆序数则也由前半段数列中 6、7 两个数来贡献，为 2。
因此 1 2 6 7 3 5 7 8 的逆序数为 4。
按这般操作，我们得到了该轮排序前数列的逆序数。而如果我们从头开始，在每一轮归并排序中都通过这样的方法来计算、累计，那么最后得到的便是原始数列的逆序数了。如果你还不太能理解整个过程，可以在纸上画画看。
我们把这个过程用代码实现，如下：
void mergeArray(int* arr, int first, int mid, int last, int* count)
{
// 新开一个动态数组，用来存已经排好的数。
// 最后会把这个动态数组存回 arr。
int* temp = (int*)malloc((last - first + 1) * sizeof(int));
int i = first, j = mid + 1, m = mid, n = last, k = 0;
while (i <= m && j <= n)
{
if (arr[i] <= arr[j])
temp[k++] = arr[i++];
else
{
temp[k++] = arr[j++];
// 贡献 (mid - i + 1) 的逆序数。
*count += mid - i + 1;
}
}
// 复制尚未到尾部的数列到结果数列。
while (i <= m)
temp[k++] = arr[i++];
while (j <= n)
temp[k++] = arr[j++];
// 存回 arr。
for (i = 0; i < k; ++i)
arr[first + i] = temp[i];
// 释放掉临时申请的动态数组。
free(temp);
}
void mergeSort(int* arr, int first, int last, int* count)
{
if (first < last)
{
int mid = (first + last) / 2;
// 对前半段数列归并排序。
mergeSort(arr, first, mid, count);
// 对后半段数列归并排序。
mergeSort(arr, mid + 1, last, count);
// 归并前后两段数列。
mergeArray(arr, first, mid, last, count);
}
}
将上述函数编写到完整程序里，可以检验其功能是否正常。
总体实现
通过上面的分析，我们已经解决了「生成全排列」和「计算逆序数」这两个关键部分。下面我们将逐步构建程序，实现整个问题的解决。
定义矩阵类型 MATRIX
为了让我们的程序拥有封装性，我们定义一种矩阵类型 MATRIX。
#define MATRIX_TYPE float
typedef struct
{
// 矩阵行数
int matHeight;
// 矩阵列数
int matWidth;
// 指向存储矩阵内容的二维指针
MATRIX_TYPE** matData;
} MATRIX;
这样，与数学中的矩阵一样，一个 MATRIX 变量即可携带行数、列数和矩阵本身的数据，这对于我们程序的迭代升级、打包封装都大有裨益。
逆序数计算函数 getReservedNumber()
利用上面我们编写的 mergeSort() 函数，我们可以写出一个封装好的，计算给定数列逆序数并返回的函数。
int getReversedNumber(int* numList, int num)
{
// 定义一个与 numList 一样长的动态数组。
int* tempList = (int*)malloc(num * sizeof(int));
// 逆序数，在函数末尾返回。
int reversedNumber = 0;
// 复制原数组到 tempList。
for (int i = 0; i < num; i++)
{
tempList[i] = numList[i];
}
// 进行归并排序，同时计算逆序数。
mergeSort(tempList, 0, num - 1, &reversedNumber);
// 释放动态数组，因为我们只需要逆序数。
free(tempList);
// 返回逆序数。
return reversedNumber;
}
在经过这一层封装之后，我们便不再需要考虑 mergeSort()、mergeArray() 等函数的相关细节了，只需要使用 getReservedNumber() 函数即可得到给定数列的逆序数，同时也不会改变给出数列的值。
完善全排列函数 makePermutation()
对我们前文所写的全排列函数 makePermutation() 进行完善。
void makePermutation(int* numList, int pHead, int pTail, MATRIX_TYPE* tempAnswer, MATRIX matMatrix)
{
if (pHead == pTail)
{
// 单次计算的答案。
MATRIX_TYPE singleAnswer = 1;
// 根据逆序数的值求取正负号。
singleAnswer *= (getReversedNumber(numList, matMatrix.matWidth) % 2 == 1) ? -1 : 1;
// 乘上 a_{1p_1}a_{2p_2}...a{np_n}。
for (int j = 0; j < matMatrix.matWidth; j++)
{
singleAnswer *= matMatrix.matData[j][numList[j]];
}
// 将单次计算的答案加到总答案上。
*tempAnswer += singleAnswer;
// 跳出函数。
return;
}
int tempExchange;
for (int i = pHead; i <= pTail; i++)
{
tempExchange = numList[pHead];
numList[pHead] = numList[i];
numList[i] = tempExchange;
makePermutation(numList, pHead + 1, pTail);
tempExchange = numList[pHead];
numList[pHead] = numList[i];
numList[i] = tempExchange;
}
}
总功能函数 matDeterminant()
将上面的一系列功能函数进行一个封装，便可以得到我们最终需要的，实现「求行列式」功能的函数——matDeterminant()。
MATRIX_TYPE matDeterminant(MATRIX matMatrix)
{
// 定义一个长度为方阵阶数的动态数组，并用 1, 2, ..., n 填充。
int* indexList = (int*)malloc(matMatrix.matHeight * sizeof(int));
for (int i = 0; i < matMatrix.matHeight; i++)
{
indexList[i] = i;
}
// 答案存储变量，从 0 开始。
MATRIX_TYPE answer = 0;
// 进行全排列。
makePermutation(indexList, 0, matMatrix.matHeight - 1, &answer, matMatrix);
// 释放掉动态数组，返回答案。
free(indexList);
return answer;
}
此外，也可以为其增加一些检测功能，来增加鲁棒性。例如，对于「求行列式」这个操作，自然要求矩阵为方阵，因此其 matHeight 和 matWidth 就应该相等。在这里省去了这些过程。