数组与特殊矩阵压缩存储

6.1 数组：线性表的推广

6.1.1 一维数组：它本来就是线性表

数组（array）是由类型相同的数据元素构成的有限序列。请注意这个定义里的两个关键词： 「类型相同」意味着每个元素占用的字节数一样，这是后面一切地址公式的前提； 「有限序列」意味着元素个数是固定的，一旦建立就不能随意增减。

一维数组 a[0..n-1] 的逻辑结构和线性表完全一致：有唯一的第一个元素、唯一的最后一个元素， 除首尾之外每个元素都有且仅有一个直接前驱和一个直接后继。所以我们可以直接说： 一维数组就是一种线性表，只是它额外约束了「元素等长」并且通常不允许插入删除。

那二维数组呢？A[m][n] 看上去是一张表格，有行有列，似乎比线性表复杂。但换个角度看： 如果我们把第 i 行的 n 个元素 a[i][0..n-1] 打包看成「一个元素」， 那么整个二维数组就退化成了一个长度为 m 的一维数组，只不过它的每个元素本身又是一个长度为 n 的线性表。换句话说：

这句话可以一直套下去：三维数组的每个元素是二维数组，二维数组的每个元素是一维数组…… 所以教材里那句「数组是线性表的推广」，说的正是这种「元素本身又是线性表」的递归结构。 理解了这一层，后面「行优先 / 列优先」不过是选择了先走哪一层而已。

6.1.2 二维数组的两种视角

图 6-1 给了我们两种看待 A[m][n] 的方式，它们直接对应两种存储方案：

行向量视角（把行当元素）

A = (R0, R1, …, Rm−1)，其中每个 Ri = (a[i][0], a[i][1], …, a[i][n−1])。 按这个顺序把元素一个个排进内存，就是行优先存储。

列向量视角（把列当元素）

A = (C0, C1, …, Cn−1)，其中每个 Cj = (a[0][j], a[1][j], …, a[m−1][j])。 按这个顺序排进内存，就是列优先存储。

两种视角在数学上是完全对称的，没有谁更「正确」。但现实中的编程语言必须选一个作为默认约定， 这就造成了 C/C++/Java/Python 与 Fortran/MATLAB 之间的差异，也造成了考试里最爱考的那一类题： 「已知某元素地址，求另一元素地址」。做题前第一件事永远是看清楚题目说的是行优先还是列优先。

6.1.3 随机存取从哪来：等长元素 + 连续存放

线性表那一章我们讲过：顺序表可以 O(1) 按下标取值，链表不行。 数组是纯粹的「顺序存储」，它的随机存取能力来自两个条件同时成立：

1. 所有元素等长：每个元素占 L 个字节，谁也不多谁也不少；
2. 所有元素连续存放：第 k 个元素紧跟在第 k−1 个元素后面，中间没有空隙。

只要这两条成立，「第 k 个元素在哪里」就不再需要「从头一个个数过去」，而可以直接乘出来：

这就是随机存取（random access）：访问任意一个元素所需的计算量完全相同， 与 k 的大小无关，所以是 O(1)（见图 6-2）。 反过来说，只要破坏其中任何一条，随机存取立刻失效： 链表的结点是 new 出来的，地址不连续，只能顺藤摸瓜， 于是按下标访问退化成 O(n)；变长字符串数组因为元素不等长，也没法用一次乘法定位。

6.1.4 数组的 ADT：为什么通常只做「取值」和「赋值」

翻开任何一本教材，数组的抽象数据类型都短得可怜，基本只有两个操作： Value(A, index) 取值、Assign(A, index, e) 赋值。 （有的书写成 GetValue / SetValue，或者干脆用 A[i] 表示。） 与之配套的还有 InitArray（初始化）、DestroyArray（销毁）， 但它们属于「生老病死」，不算常规操作。

为什么数组不像线性表那样提供 Insert / Delete？原因有三条，一条比一条本质：

• 结构一旦建立就固定不变。数组的维数和各维长度在初始化时就确定了， A[m][n] 的含义是「恰好 m×n 个元素」。插入一个元素会让它变成 m×n+1 个，这个形状根本不是一个合法的二维数组。
• 插入会摧毁随机存取。如果允许在中间插一个元素、后面的元素整体后移， 那么每次插入都是 O(n)；更糟的是，一旦允许元素个数变化， 「连续 + 等长」的假设就要靠不断搬迁来维持，数组最大的优势就没了。
• 数组的定位是「底层存储结构」。数组不是给人直接用的高层容器， 而是实现别的结构的砖头。顺序表要在中间插入时，会用数组 + 元素个数计数 + 容量 自己管理：把「有多少个有效元素」这件事交给一个额外的 length 变量， 而不是交给数组。这样数组本身保持纯粹，随机存取能力得以保留。

6.2 多维数组的存储：行优先与列优先

6.2.1 行优先（Row-Major）：C/C++ 的选择

二维数组在逻辑上是二维的，但计算机内存是一维的，一根长长的地址轴。 所谓「存储」，就是把这 m×n 个元素排成一个队列塞进这根地址轴里。 排列方式不止一种，最常见的就是行优先（row-major order）： 先存第 0 行，再存第 1 行……同一行里从左到右依次存。

现在推导地址公式。设 A[m][n]，首元素 a[0][0] 的地址是 LOC(a[0][0])，每个元素占 L 个字节。要求 a[i][j] 的地址， 只需要回答一个问题：排在 a[i][j] 前面的元素一共有多少个？

1. 第 0 行到第 i−1 行都是完整地排在它前面的，一共 i 行， 每行 n 个元素，所以有 i × n 个；
2. 在它自己所在的第 i 行里，列号比 j 小的还有 j 个元素 （列号从 0 到 j−1）；
3. 合计前面有 i × n + j 个元素，每个占 L 字节， 所以地址要在首地址上加上 (i × n + j) × L。

公式里的 i × n + j 其实就是元素 a[i][j] 在「拉直成一行」之后的下标， 我们把它叫做线性下标或一维序号，记作 k = i × n + j。 考试时如果你一时想不起公式，最快的办法就是举个小例子现场数： 比如 A[3][4] 里的 a[1][2]，前面有 1 整行（4 个）+ 本行 2 个 = 6 个， 于是 k = 6，与 1×4+2 = 6 完全一致。

6.2.2 列优先（Column-Major）：Fortran / MATLAB 的选择

列优先的排法刚好反过来：先存第 0 列，再存第 1 列，同一列里从上到下依次存。 推导完全对称，只是「整块」的单位从「行」变成了「列」：

1. 第 0 列到第 j−1 列都是完整的 j 列，每列 m 个元素，共 j × m 个；
2. 本列中行号比 i 小的还有 i 个；
3. 合计 j × m + i 个元素排在前头。

请特别留意：行优先的公式里出现的是列数 n，列优先的公式里出现的是行数 m。 这是最容易写错的地方，因为直觉上我们总觉得「求 (i,j) 的地址应该用到 i 前面的东西、乘上列数」， 但在列优先里，跨一整列要跨过的是 m 个元素。记忆口诀：行优先「乘列数」，列优先「乘行数」。

下面这段代码把两种排布都打印出来，可以直接对照上面两张图看：

#include <iostream>
using namespace std;

/* ============================================================
   二维数组的地址计算：行优先 vs 列优先
   设 A 是 m 行 n 列的二维数组，第一个元素 a[0][0] 的地址是 base，
   每个元素占 L 个字节（所有元素等长）。

   行优先 Row-Major（C / C++ / Java / Python 采用）：
       LOC(a[i][j]) = base + (i * n + j) * L
   列优先 Column-Major（Fortran / MATLAB 采用）：
       LOC(a[i][j]) = base + (j * m + i) * L
   ============================================================ */

long long rowMajorAddr(long long base, int n, int i, int j, int L) {
    return base + (1LL * i * n + j) * L;      /* 前面有 i 整行，本行再走 j 格 */
}

long long colMajorAddr(long long base, int m, int i, int j, int L) {
    return base + (1LL * j * m + i) * L;      /* 前面有 j 整列，本列再走 i 格 */
}

int main() {
    const int m = 3, n = 4, L = 4;            /* A[3][4]，每个 int 占 4 字节 */
    const long long base = 1000;              /* 假设 LOC(a[0][0]) = 1000 */

    cout << "A[3][4], base = " << base << ", L = " << L << "\n";

    cout << "row-major order: ";
    for (int i = 0; i < m; ++i)
        for (int j = 0; j < n; ++j)
            cout << "a[" << i << "][" << j << "]@" << rowMajorAddr(base, n, i, j, L) << "  ";
    cout << "\n";

    cout << "col-major order: ";
    for (int j = 0; j < n; ++j)
        for (int i = 0; i < m; ++i)
            cout << "a[" << i << "][" << j << "]@" << colMajorAddr(base, m, i, j, L) << "  ";
    cout << "\n";

    /* 单独验算 a[1][2]：行优先排在它前面的元素是 1*4+2 = 6 个 */
    cout << "a[1][2] row-major: (1*4+2)*4 = " << (1 * n + 2) * L
         << " bytes, addr = " << rowMajorAddr(base, n, 1, 2, L) << "\n";
    /* 列优先排在它前面的是 2*3+1 = 7 个 */
    cout << "a[1][2] col-major: (2*3+1)*4 = " << (2 * m + 1) * L
         << " bytes, addr = " << colMajorAddr(base, m, 1, 2, L) << "\n";

    /* 经典考题：题目不给首地址，只给某个元素的地址，问另一个元素的地址。
       做法是"只算偏移量之差"，base 会在相减时消掉。 */
    long long loc32 = rowMajorAddr(base, n, 3, 2, L);              /* 已知 a[3][2] 的地址 */
    long long delta = (1LL * 1 * n + 3) - (1LL * 3 * n + 2);       /* a[1][3] 与 a[3][2] 的元素序号差 */
    cout << "given a[3][2]@" << loc32 << ", a[1][3]@" << loc32 + delta * L
         << " (delta = " << delta << " elements)\n";
    return 0;
}

6.2.3 推广到 k 维数组：从最后一维开始递推

现在把二维的结论推广到 k 维。设数组是 d1 × d2 × … × dk， 下标记作 (i1, i2, …, ik)，各维下标都从 0 开始。 在行优先（最后一维变化最快）的约定下，线性序号是：

这个公式不要死记，用递推（从最后一维往前套）的方式想，三秒钟就能重建出来：

• 最后一维 ik 每加 1，就在内存里往后走 1 个元素，所以它的权重是 1；
• 倒数第二维 ik−1 每加 1，意味着「把最后一维整段走完一遍」， 所以它的权重是 dk；
• 再往前一维 ik−2 每加 1，意味着「把后面两维整块走完一遍」， 权重是 dk−1 × dk；
• 依此类推，第 t 维的权重就是「它后面所有维度长度之积」。

用二维验证一下：d1=m, d2=n，则 i1 的权重是 d2 = n， i2 的权重是 1，于是 k = i×n + j，正是行优先公式，完全吻合。 如果换成列优先，规律反过来：第一维变化最快，所以第一维权重为 1，越靠前的维度权重越小。

#include <iostream>
#include <vector>
using namespace std;

/* ============================================================
   k 维数组的行优先地址
   设数组是 d1 x d2 x ... x dk，下标 (i1, i2, ..., ik) 均从 0 开始，
   则线性序号（第几个元素）= i1*(d2*d3*...*dk) + i2*(d3*...*dk) + ... + ik
   LOC = base + 线性序号 * L

   递推理解：最后一维 ik 每走一步只跨 1 个元素；
             倒数第二维 i(k-1) 每走一步要跨 dk 个元素；
             再往前一层，每走一步跨 dk*d(k-1) 个……
             于是从最后一维开始往前"套"，就得到了上面的公式。
   ============================================================ */

long long kdLinearIndex(const vector<int>& dim, const vector<int>& idx) {
    long long lin = 0;
    for (size_t t = 0; t < dim.size(); ++t) {
        long long stride = 1;                       /* 第 t 维的步长 = 它后面所有维度的乘积 */
        for (size_t u = t + 1; u < dim.size(); ++u) stride *= dim[u];
        lin += 1LL * idx[t] * stride;
    }
    return lin;
}

/* 反向操作：已知线性序号，拆回各维下标（从最后一维往前取余、取商） */
vector<int> kdDecode(const vector<int>& dim, long long lin) {
    vector<int> idx(dim.size(), 0);
    for (size_t t = dim.size(); t-- > 0; ) {
        idx[t] = (int)(lin % dim[t]);
        lin /= dim[t];
    }
    return idx;
}

int main() {
    vector<int> dim = {3, 4, 5};                 /* 三维数组 A[3][4][5] */
    vector<int> idx = {1, 2, 3};                 /* 求 a[1][2][3] */
    const long long base = 1000, L = 4;

    long long lin = kdLinearIndex(dim, idx);
    cout << "A[3][4][5] a[1][2][3]: linear index = "
         << "1*(4*5) + 2*5 + 3 = " << lin << "\n";                 /* 20 + 10 + 3 = 33 */
    cout << "  LOC = " << base << " + " << lin << "*" << L << " = "
         << base + lin * L << "\n";                                /* 1132 */

    vector<int> back = kdDecode(dim, lin);
    cout << "decode(" << lin << ") = [" << back[0] << ", " << back[1]
         << ", " << back[2] << "]\n";

    /* 列优先（第一维变得最快）：线性序号 = i1 + i2*d1 + i3*d1*d2 */
    long long colLin = idx[0] + 1LL * idx[1] * dim[0] + 1LL * idx[2] * dim[0] * dim[1];
    cout << "a[1][2][3] col-major: 1 + 2*3 + 3*3*4 = " << colLin
         << ", LOC = " << base + colLin * L << "\n";               /* 43, 1172 */

    /* 相邻元素验证：把整块内存按行优先铺开，看看每个位置属于哪个下标 */
    cout << "row-major memory layout of A[2][3][2]:\n";
    vector<int> d2 = {2, 3, 2};
    for (long long p = 0; p < 2 * 3 * 2; ++p) {
        vector<int> ij = kdDecode(d2, p);
        cout << p << ":a[" << ij[0] << "][" << ij[1] << "][" << ij[2] << "]  ";
    }
    cout << "\n";
    return 0;
}

6.2.4 四道地址计算例题（务必自己先算一遍）

例 6-1 设二维数组 A[4][5]（4 行 5 列），每个元素占 4 个字节， LOC(a[0][0]) = 1000。求 a[2][3] 在行优先与列优先下的地址。

例 6-2 设 A[5][6] 按行优先存储，每个元素 2 字节， LOC(a[0][0]) = 2000。求 a[3][4] 与 a[4][5] 的地址， 并问：从 a[3][4] 到 a[4][5] 跨过了多少个元素？

例 6-3 设三维数组 A[3][4][5]，每个元素 4 字节，LOC(a[0][0][0]) = 1000。 求 a[1][2][3] 在行优先下的地址；若改为列优先，地址又是多少？

例 6-4 设二维数组 A[10][20] 按行优先存储，每个元素占 1 个字节。 已知 a[3][5] 的地址是 1010，求 a[5][8] 的地址。

6.2.5 动画演示：同一份数据，两种压平方式

下面这个动画把 A[3][4] 的 12 个格子逐个「搬」到底下的内存带上。 先按行优先走一遍，再按列优先走一遍，请留意每一步里 二维下标 → 线性下标 k → 内存偏移 k×L 这三者的对应关系：

6.3 矩阵的压缩存储：压什么、怎么找回来

6.3.1 什么叫压缩存储

矩阵是数值计算的主角，一个 m × n 的矩阵用二维数组存下来要占 m × n × L 个字节。1000×1000 的 double 矩阵就是 8 MB， 听起来还行；但工程里动辄几十万阶的矩阵，直接开二维数组会瞬间把内存吃光。 于是我们自然会问：这些空间里，有多少是在存「废话」？

答案是很多。看几个例子：

• 对称矩阵：a[i][j] = a[j][i]，以主对角线为轴的两半完全相同， 存两遍等于白白浪费一半空间。
• 三角矩阵：主对角线的一侧全是同一个常数 c（常见于 c = 0）， 这一侧成千上万个格子存的都是同一个数。
• 三对角（带状）矩阵：除了三条对角线，其余全是 0。n 阶矩阵有 n² 个格子， 真正有用的只有 3n−2 个。
• 稀疏矩阵：非零元素极少（通常不到 5%），绝大多数格子是 0。

针对这些情况，教材给出了压缩存储的定义：

这个定义里有两个关键词值得抠一抠。第一是「多个值相同的元素只分配一个空间」—— 注意它说的是「值相同」，而不是「位置相邻」。对称矩阵里 a[i][j] 与 a[j][i] 位置上离得很远，但值一定相同，于是我们可以只留一个。 第二是「零元素不分配空间」——注意这在数学上是有前提的： 只有当「0 是这批数据的默认值」时，不存 0 才能靠「没存的就是 0」把信息补回来。 如果矩阵里 0 和别的数一样是「有意义的数据」，那就不能这么干。

6.3.2 压缩的前提：元素分布必须有规律

压缩存储不是变魔术，它必须付出代价。代价是什么？是「要知道某个位置的元素在哪里， 得先算一算」。而这个「算」能不能在 O(1) 时间内完成， 完全取决于元素分布的规律性。

所以「压缩存储」实际上分成两大流派： 公式法（顺序存储 + 下标映射）用于对称 / 三角 / 带状这类结构规整的矩阵； 列表法（三元组 / 链表）用于稀疏矩阵。前者能保住 O(1) 随机存取， 后者不能。这也正是本章 6.4–6.6 与 6.7 的分界线。

6.3.3 下标映射函数：压缩存储的灵魂

压缩之后，原来的二维下标 (i, j) 不再等于内存位置，我们需要一个函数把它翻译过去：

这个 f 就叫下标映射函数（index mapping function）。有了它， get(i,j) 就是 sa[f(i,j)]，set(i,j,v) 就是 sa[f(i,j)] = v，两条语句、常数时间。整个压缩存储方案的价值， 全都压在这个函数身上。

好消息是：所有这类函数的推导套路都是同一套三步法，掌握之后可以临场推、不用背：

1. 第一步：确定「存哪些」。画出矩阵，把要存的区域涂上颜色， 用下标不等式描述它（比如对称矩阵存 i ≥ j）。 顺便数一数总共有多少个元素，得到数组长度。
2. 第二步：确定「按什么顺序存」。绝大多数教材用行优先： 先存第 0 行里要存的那些，再存第 1 行……这一步决定了「整块」是行还是列。
3. 第三步：数「排在 (i,j) 前面的元素有多少个」。 把它拆成「前面完整的若干行」+「本行内排在它前面的几个」， 前者是一个等差数列求和，后者往往就是几。加在一起就是 k。

6.3.4 压完之后还能随机存取吗

这是考试里最爱问的辨析题，答案要看「f 能不能 O(1) 算出来」：

一句话总结：压缩存储是否保留随机存取，取决于压缩后「位置信息」是否被一个 O(1) 可算的公式完整保留下来。对称 / 三角 / 带状矩阵的规律是「全局的」， 一个公式管全部；稀疏矩阵的规律是「局部的」，只能靠额外的索引结构（rpos[]、 十字链表、哈希）去补救。

6.4 对称矩阵的压缩存储

6.4.1 定义与压缩方案

设 A 是 n 阶方阵，如果对所有 0 ≤ i, j ≤ n−1 都有

则称 A 为对称矩阵（symmetric matrix）。它是沿主对角线「照镜子」左右对称的， 所以只要知道一半，另一半就全知道了。

压缩方案的选择很自然：只存下三角（含主对角线），按行优先顺序 压进一个一维数组 sa[]。为什么不存上三角？纯粹是习惯—— 教材和考试统一按下三角来，你跟着来就行（存上三角也完全可以，公式对称地改一下即可）。

为什么是「含主对角线」？因为对角线上的元素 a[i][i] 满足 i = j， 它在镜子上的像是它自己，只存一次正好，丢掉就再也找不回来了。

6.4.2 下标映射公式的完整推导

现在按 6.3.3 的三步法，把 k = f(i, j) 推出来。这一段请务必跟着算一遍， 它是本章所有公式的模板。

第一步：存哪些？共多少个？

存下三角（含对角线），即所有满足 i ≥ j 的元素。数一数每行有几个： 第 0 行有 1 个（a[0][0]），第 1 行有 2 个（a[1][0], a[1][1]）， …，第 i 行有 i+1 个，…，第 n−1 行有 n 个。 总数是一个等差数列之和：

用等差数列求和公式（首项 + 末项）× 项数 ÷ 2：

所以一维数组 sa[] 的长度是 n(n+1)/2， 下标范围是 0 … n(n+1)/2 − 1。 例如 n = 5 时长度是 5×6/2 = 15，与图 6-5 一致。 （顺带记一个常用结论：1+2+…+i = i(i+1)/2， 这个「三角形数」在后面反复出现，是本章最常用的求和式。）

第二步：按什么顺序存？

按行优先——先把第 0 行里要存的元素从左到右放好，再放第 1 行的，依此类推。 于是在一维数组里，第 0 行占据 sa[0..0]， 第 1 行占据 sa[1..2]，第 2 行占据 sa[3..5]， 第 i 行占据 sa[i(i+1)/2 .. i(i+1)/2 + i]。 这正是图 6-5 右边那个「阶梯」形状的来历。

第三步：数一数 (i, j) 前面有多少个元素。

分两种情况讨论。

情况 A：i ≥ j（元素本来就在下三角里，值就存在 sa 中）。 排在它前面的元素由两部分组成：

• 第 0 行到第 i−1 行是完整的 i 行， 每行的元素个数分别是 1, 2, …, i， 共 1 + 2 + … + i = i(i+1)/2 个；
• 在它自己所在的第 i 行里，列号从 0 数到 j−1， 共 j 个元素（注意：是 j 个，不是 j+1 个， 因为 a[i][j] 自己不算「排在前面」）。

两部分相加：

情况 B：i < j（元素在上三角，本身不存储）。 由对称性 a[i][j] = a[j][i]，而 (j, i) 满足 j > i， 即 j ≥ i，一定落在下三角里（含对角线），可以直接套用情况 A 的公式， 只要把 i 和 j 的角色对调：

把两个式子合起来写，就是对称矩阵的下标映射函数：

6.4.3 两道手算例题

例 6-5 n = 5 的对称矩阵按「下三角 + 行优先」压缩进一维数组 sa[]， 问 a[3][1] 与 a[1][3] 分别存在 sa 的哪个位置？ 数组总长是多少？

例 6-6 n = 6 的对称矩阵压缩后，a[5][5]、a[2][5]、 a[4][2] 分别存在哪里？一维数组最大下标是多少？

例 6-7（反向题） 一个 6 阶对称矩阵压缩后得到 sa[0..20]。已知某个值存在 sa[12]， 它在原矩阵中的下标 (i, j) 是什么（只写下三角的那个）？

6.4.4 C++ 实现：全局数组 + 下标映射函数

把公式写进代码只有两行的事。竞赛里的写法非常直接：开一个全局数组 sa[]，再写一个 idx(i, j) 把二维下标映射成一维下标， 存取各写一个自由函数就行——不封类、不模板、不抛异常，代码能在考场上三分钟敲完。 要注意的只有两点：下标越界的检查（越界只提示、不崩溃）和 映射函数要单独写出来，方便和手算的 k 对照。

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   对称矩阵的压缩存储（竞赛写法：全局数组 + 下标映射函数）
   定义：n 阶方阵 A 满足 a[i][j] == a[j][i]，只需存下三角（含对角线）。

   按行优先把下三角压进一维数组 sa[0 .. n(n+1)/2 - 1]：
       i >= j （下三角）：k = i*(i+1)/2 + j
       i <  j （上三角）：k = j*(j+1)/2 + i     （利用 a[i][j] = a[j][i]）

   存储长度 1 + 2 + ... + n = n(n+1)/2，比原来的 n^2 省了近一半。
   例如 n = 5：长度 15，a[3][1] 与 a[1][3] 都映射到 k = 7，共用同一个单元。
   ============================================================ */

const int N  = 1005;                  /* 最大阶数 */
const int SA = N * (N + 1) / 2;       /* 一维数组长度：n(n+1)/2 个单元 */
int sa[SA];                           /* 压缩数组：行优先存下三角（含对角线） */
int n;                                /* 矩阵阶数 */

/* 下标映射函数：把二维下标 (i, j) 映射成一维下标 k
   下三角：前面有 i 整行共 1+2+...+i = i(i+1)/2 个元素，本行再走 j 格
   上三角：换成对称位置 (j, i) 套同一个公式 —— 「对称性免费」就是这么来的 */
inline int idx(int i, int j) {
    return (i >= j) ? (i * (i + 1) / 2 + j)      /* 下三角 */
                    : (j * (j + 1) / 2 + i);     /* 上三角 */
}

/* 边界检查：竞赛里不抛异常，越界就打印提示并拒绝这次操作 */
bool inRange(int i, int j) {
    if (0 <= i && i < n && 0 <= j && j < n) return true;
    printf("  [warn] (%d, %d) out of range\n", i, j);
    return false;
}

/* 存：只需写一次，对称位置「自动」也是同一个值 */
void setVal(int i, int j, int v) {
    if (!inRange(i, j)) return;
    sa[idx(i, j)] = v;
}

/* 取：上、下三角走同一套映射，O(1)（一次乘法 + 一次加法） */
int getVal(int i, int j) {
    if (!inRange(i, j)) return 0;
    return sa[idx(i, j)];
}

/* 打印整个一维数组，方便与手算结果对照 */
void dump() {
    int len = n * (n + 1) / 2;
    for (int k = 0; k < len; ++k) printf("%d%c", sa[k], k + 1 == len ? '\n' : ' ');
}

/* 函数名为什么叫 setVal 而不是 set？
   using namespace std 之后，自由函数 set 会和 std::set 撞名，
   g++ 直接报 "reference to 'set' is ambiguous"，所以换个名字。 */

int main() {
    n = 5;
    printf("n = %d, storage size = %d (n(n+1)/2 = %d)\n",
           n, n * (n + 1) / 2, n * (n + 1) / 2);          /* 15 */

    /* 手算例题：n = 5 时 a[3][1] 与 a[1][3] 存在哪个位置？ */
    printf("a[3][1] -> k = %d\n", idx(3, 1));             /* 3*4/2+1 = 7 */
    printf("a[1][3] -> k = %d\n", idx(1, 3));             /* 3*4/2+1 = 7 */

    /* 按下三角的行优先顺序填入 1,2,3,...，再打印一维数组对照 */
    int val = 0;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j <= i; ++j)
            setVal(i, j, ++val);
    printf("sa[] = ");
    dump();                                               /* 1 2 3 ... 15 */

    /* 验证对称性：任意 (i,j) 都有 getVal(i,j) == getVal(j,i) */
    bool ok = true;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            if (getVal(i, j) != getVal(j, i)) ok = false;
    printf("symmetry check: %s\n", ok ? "PASS" : "FAIL");

    /* 写一次 a[3][1]，读 a[1][3] 也是同一个值：两次访问命中同一个内存单元 */
    setVal(3, 1, 41);
    printf("set a[3][1] = 41, then a[1][3] = %d\n", getVal(1, 3));   /* 41 */

    /* 边界检查：越界只提示、不崩溃 */
    printf("a[5][0] = %d  (out of range)\n", getVal(5, 0));
    return 0;
}

这段代码里有三个细节值得单独说：

• idx(i, j) 用 (i >= j) ? ... : ... 一次搞定上下三角， 这就是公式的直接翻译。它只需要一次乘法和一次加法，所以 getVal/setVal 是 O(1)——对称矩阵压缩后随机存取能力完好无损。
• setVal(3,1,41) 之后 getVal(1,3) 也是 41，这不是巧合， 而是因为 idx(3,1) 与 idx(1,3) 都等于 7， 两次访问命中同一个内存单元。「对称性」在压缩存储里是免费的。
• 数组开成全局的 sa[N*(N+1)/2]：竞赛里数组开够就不再动态申请内存， 程序结束由操作系统统一回收，不写 new/delete，也不写析构函数。 代价是长度上限写死在 N 里，做题时按题目数据范围开即可。
• 写函数时别取名 set：using namespace std 之后它会和 std::set 撞名，g++ 直接报 reference to 'set' is ambiguous， 所以这里叫 setVal（同理 getVal，避免和 std::get 混淆）。

6.4.5 动画演示：映射过程与对称性

下面的动画把 5 阶对称矩阵的下三角逐个搬进 sa[]。 每次处理 (i, j) 时，画面会同时标出它的「镜像位置」(j, i)， 并显示两者算出来的 k 完全相同——这就是对称矩阵能省一半空间的直观解释：

6.5 三角矩阵的压缩存储

6.5.1 两种三角矩阵与压缩方案

三角矩阵（triangular matrix）分两种，定义只差一个方向：

下三角矩阵（lower triangular）

主对角线以上（不含主对角线）的元素全是同一个常数 c， 即当 i < j 时恒有 a[i][j] = c。 主对角线及以下的部分是「真实数据」。

上三角矩阵（upper triangular）

主对角线以下（不含主对角线）的元素全是同一个常数 c， 即当 i > j 时恒有 a[i][j] = c。

注意定义里那句「不含主对角线」——对角线永远是真实数据，必须老老实实存下来。 最常见的常数是 c = 0（例如高斯消元过程中出现的上三角方程组系数矩阵）， 但定义本身对 c 没有限制，所以我们的方案必须能存下任意常数。

压缩方案与对称矩阵几乎一样（两种矩阵的对照见图 6-6），只有一处不同：

1. 把真实数据区（下三角或上三角，含对角线）按行优先压进 sa[0 … n(n+1)/2 − 1]，共 n(n+1)/2 个单元；
2. 常数 c 单独占用最后一个单元 sa[n(n+1)/2]。

所以三角矩阵的压缩数组长度是 n(n+1)/2 + 1。 为什么常数区已经「全都一样」了，还要单独留一个位置？因为「全都一样」不等于 「等于多少已知」。对称矩阵的上三角可以靠 a[j][i] 推算出来， 而三角矩阵的常数区推算不出来，必须把 c 的值存一份。 （如果题目明确说了 c = 0，实现时当然可以省掉这个单元， 但考试里请按通式 +1 回答。）

6.5.2 下三角矩阵：公式与长度

下三角矩阵要存的部分正好是 i ≥ j，和对称矩阵存的区域一模一样， 所以下标公式也一模一样：

唯一的区别出在 i < j 的时候：对称矩阵会「照镜子」去查 a[j][i]， 而下三角矩阵里这些位置的值是常数 c，直接返回 sa[n(n+1)/2] 即可。 写成代码就是一句话：

数组长度 n(n+1)/2 + 1：前 n(n+1)/2 个单元放真实数据， 最后 1 个单元放常数。

6.5.3 上三角矩阵：公式推导

上三角矩阵的推导要用到 6.3.3 的同一个套路，只是「每行有几个元素」不再是 1,2,3,…， 而是递减的：

第一步：存哪些、共多少个？

按行优先存 i ≤ j 的部分。第 i 行要存的是 j = i, i+1, …, n−1，个数是 n − i。 于是各行个数依次是 n, n−1, n−2, …, 1，总数仍是

第二步：按行优先存。第 i 行的起点是「前 i 行的总数」。

第三步：数 (i, j) 前面有多少个元素。分成两块：

• 前 i 行的总数：n + (n−1) + … + (n−i+1)。 这是一个首项为 n、末项为 n−i+1、共 i 项的等差数列， 用「(首项 + 末项) × 项数 ÷ 2」：
  (n + (n−i+1)) × i / 2 = i(2n − i + 1)/2
  （也可以用另一种算法验证：Σt=0i−1(n−t) = i·n − i(i−1)/2， 通分后是 (2in − i² + i)/2 = i(2n−i+1)/2，两式一致 ✅）
• 本行内排在 j 前面的个数：本行从列号 i 开始放， 放到 j 之前一共放了 j − i 个。

两式相加，得到上三角矩阵的下标映射函数：

当 i > j 时，元素落在常数区，k = n(n+1)/2。 数组长度同样是 n(n+1)/2 + 1。

6.5.4 两道手算例题

例 6-8 n = 5 的下三角矩阵（常数 c = 0）按行优先压缩， 求 a[3][2]、a[4][4]、a[0][4] 的存储位置， 以及数组总长度。

例 6-9 同样 n = 5，但换成上三角矩阵（常数 c = 0）。 求 a[1][3]、a[0][4]、a[4][4]、a[3][1] 的位置。

6.5.5 C++ 实现：两套映射函数与边界检查

下三角与上三角各给一套映射函数（lowIdx / upIdx）和两个数组 （saL[] / saU[]）。除了公式本身，请重点看边界检查是怎么做的： 下标越界要提示；往常数区写一个不等于 c 的值，也要立刻报错（这里用返回 false 表示拒绝）——因为那说明调用者的逻辑已经错了，静默接受只会让 bug 藏得更深。 竞赛代码不抛异常，用返回值或打印提示来表达「这次操作不合法」，一样能起到暴露错误的作用。

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   三角矩阵的压缩存储（竞赛写法：两套映射函数 + 常数单元单独存）
   ------------------------------------------------------------
   下三角矩阵：主对角线以上（不含对角线）的元素全是同一个常数 c
       存法：saL[0 .. n(n+1)/2 - 1] 按下三角行优先存；
             saL[n(n+1)/2] 单独存常数 c。
       i >= j : k = i*(i+1)/2 + j
       i <  j : k = n(n+1)/2             （常数单元）
       数组长度 n(n+1)/2 + 1

   上三角矩阵：主对角线以下（不含对角线）的元素全是同一个常数 c
       存法：按行优先存上三角，第 i 行存 j = i .. n-1，共 n-i 个元素；
             前 i 行共 sum_{t=0}^{i-1}(n-t) = i*n - i(i-1)/2 = i(2n-i+1)/2 个元素。
             注意第 i 行是从列 i 开始存的，所以行内偏移是 (j-i) 而不是 j。
       i <= j : k = i*(2n-i+1)/2 + (j-i)
       i >  j : k = n(n+1)/2
       数组长度同样 n(n+1)/2 + 1
   ============================================================ */

const int N  = 1005;                  /* 最大阶数 */
const int SA = N * (N + 1) / 2 + 1;   /* 长度：n(n+1)/2 个真实数据 + 1 个常数单元 */
int saL[SA];                          /* 下三角矩阵的压缩数组 */
int saU[SA];                          /* 上三角矩阵的压缩数组 */
int n;                                /* 阶数 */
int cst;                              /* 常数 c */

/* 常数单元的下标：下三角、上三角都是 n(n+1)/2 这一格 */
inline int constIdx() { return n * (n + 1) / 2; }

/* 下三角映射（i >= j 时有效）：k = i*(i+1)/2 + j */
inline int lowIdx(int i, int j) { return i * (i + 1) / 2 + j; }

/* 上三角映射（i <= j 时有效）：k = i*(2n-i+1)/2 + (j-i) */
inline int upIdx(int i, int j) { return i * (2 * n - i + 1) / 2 + (j - i); }

/* 边界检查：竞赛里不抛异常，越界就返回 false 让调用者自己处理 */
bool inRange(int i, int j) {
    if (0 <= i && i < n && 0 <= j && j < n) return true;
    printf("  [warn] (%d, %d) out of range\n", i, j);
    return false;
}

/* 初始化：把两个数组的常数单元都预置成 c，
   这样即使从没写过常数区，读出来的也是正确的 c */
void init(int nn, int c) {
    n = nn; cst = c;
    for (int k = 0; k <= constIdx(); ++k) { saL[k] = c; saU[k] = c; }
}

/* ---------- 下三角 ---------- */
/* 存：i >= j 写真实位置；i < j 落在常数区，只允许写 c */
bool lowSet(int i, int j, int v) {
    if (!inRange(i, j)) return false;
    if (i < j) {                     /* 这一类位置压缩后根本没有独立存储单元 */
        if (v != cst) {
            printf("  [warn] a[%d][%d] 在常数区，只能写 c = %d\n", i, j, cst);
            return false;            /* 静默接受只会让 bug 藏得更深 */
        }
        return true;                 /* 写的就是 c，等于没写 */
    }
    saL[lowIdx(i, j)] = v;
    return true;
}
int lowGet(int i, int j) {
    if (!inRange(i, j)) return 0;
    return i >= j ? saL[lowIdx(i, j)] : saL[constIdx()];
}

/* ---------- 上三角 ---------- */
bool upSet(int i, int j, int v) {
    if (!inRange(i, j)) return false;
    if (i > j) {                     /* 常数区 */
        if (v != cst) {
            printf("  [warn] a[%d][%d] 在常数区，只能写 c = %d\n", i, j, cst);
            return false;
        }
        return true;
    }
    saU[upIdx(i, j)] = v;
    return true;
}
int upGet(int i, int j) {
    if (!inRange(i, j)) return 0;
    return i <= j ? saU[upIdx(i, j)] : saU[constIdx()];
}

int main() {
    init(5, 0);                                            /* n = 5，常数 c = 0 */

    lowSet(3, 2, 7);
    lowSet(4, 4, 9);
    printf("lower triangular, n = %d\n", n);
    printf("  a[3][2] = %d, k = %d   (3*4/2 + 2 = 8)\n", lowGet(3, 2), lowIdx(3, 2));
    printf("  a[4][4] = %d, k = %d  (4*5/2 + 4 = 14)\n", lowGet(4, 4), lowIdx(4, 4));
    printf("  a[0][4] = %d, k = %d  (constant cell)\n", lowGet(0, 4), constIdx());
    printf("  storage size = %d = n(n+1)/2 + 1 = %d\n",
           constIdx() + 1, n * (n + 1) / 2 + 1);           /* 16 */

    upSet(1, 3, 5);
    upSet(0, 4, 6);
    printf("upper triangular, n = %d\n", n);
    printf("  a[1][3] = %d, k = %d   (1*(10-1+1)/2 + 2 = 7)\n", upGet(1, 3), upIdx(1, 3));
    printf("  a[0][4] = %d, k = %d   (0*11/2 + 4 = 4)\n", upGet(0, 4), upIdx(0, 4));
    printf("  a[4][4] = %d, k = %d  (4*(10-4+1)/2 + 0 = 14)\n", upGet(4, 4), upIdx(4, 4));
    printf("  a[4][0] = %d  (constant cell)\n", upGet(4, 0));
    printf("  storage size = %d\n", constIdx() + 1);

    /* 边界检查：越界、以及「往常数区写非 c 值」都应该被拦住 */
    printf("lowGet(5, 0)      -> %d\n", lowGet(5, 0));
    printf("lowSet(0, 4, 1)   -> %s\n", lowSet(0, 4, 1) ? "true" : "false");
    return 0;
}

三个要点：

• inRange(i, j) 负责越界检查。竞赛里不用 throw std::out_of_range 那套异常机制（考场上没人写 try/catch），而是打印一行提示并拒绝这次操作： 越界属于调用者的编程错误，必须立刻看得见；打印出来的 warning 就是最直接的暴露方式。
• lowSet(i, j, v) 在 i < j（下三角的常数区）时， 如果 v != cst 直接返回 false： 这一类位置在压缩后根本没有独立的存储单元，允许写入不同值会让数据结构自相矛盾。
• init(n, c) 把两个数组的常数单元都预置成 c （位置就是 saL[n(n+1)/2] / saU[n(n+1)/2]）， 这样即使从没写过常数区，读出来的也是正确的 c—— 这就是「三角矩阵要比对称矩阵多存一个单元」的那一格。

6.6 对角矩阵（带状矩阵）的压缩存储

6.6.1 三对角矩阵的定义

如果矩阵中所有非零元素都集中在主对角线附近的若干条对角线上，就称它为 带状矩阵（band matrix），集中分布的这些对角线构成「带」。 最常见的是三对角矩阵（tridiagonal matrix），也叫带宽为 3 的带状矩阵：

三条对角线各有名字：i = j 是主对角线，j = i + 1 是主对角线 上方的第一条，i = j + 1 是主对角线下方的第一条。 三对角矩阵在数值计算里出现的频率高得惊人：三次样条插值、一维热传导方程的差分格式、 隐式求解偏微分方程，最后都会归结为「解一个三对角线性方程组」， 而这类方程组可以用追赶法（Thomas 算法）在 O(n) 时间内解出—— 前提正是它被压缩存储了。

一个 n 阶三对角矩阵有 n² 个元素，但真正可能非零的只有多少？ 数一数：第 0 行 2 个；中间第 1 … n−2 行每行 3 个；最后一行 2 个。于是

例如 n = 1000 时，n² = 1000000 而 3n − 2 = 2998， 压缩率高达 99.7%。这就是压缩存储的威力。

6.6.2 下标公式的推导与长度 3n−2

还是那三步。这次「每行存几个」不再是等差的了，所以要分情况。

第一步：每行存几个？

• 第 0 行：只有 a[0][0] 和 a[0][1]，共 2 个 （a[0][−1] 不存在）；
• 第 i 行（1 ≤ i ≤ n−2）：a[i][i−1], a[i][i], a[i][i+1]，共 3 个；
• 第 n−1 行：a[n−1][n−2] 和 a[n−1][n−1]，共 2 个。

总数：

于是压缩数组 sa[] 的长度是 3n − 2，下标范围 0 … 3n−3。 n = 5 时是 13，n = 6 时是 16（与图 6-7 一致）。

第二步与第三步：算 k。

情况 1：i = 0。第 0 行前面什么都没有，本行从列 0 开始放， 所以行内偏移就是 j：

情况 2：i ≥ 1。前 i 行（第 0 行到第 i−1 行）的元素总数是 「第 0 行的 2 个 + 中间 (i−1) 行每行 3 个」：

第 i 行是从列 i−1 开始放的（因为第一个元素是 a[i][i−1]）， 所以 a[i][j] 在本行内是第 j − (i−1) = j − i + 1 个元素， 行内偏移就是 j − i + 1。两者相加：

有意思的是，i = 0 时 2i + j = j，与情况 1 的结果完全一致。 所以两种写法可以统一成一个式子：

为什么能统一？因为第 0 行「少的那一个元素」正好被 3i − 1 在 i = 0 时 多减的 1 抵消掉了。这也解释了一个容易困惑的点：这个统一公式只在带内成立， 带外的 (i, j) 代进去会算出越界的 k（比如 a[0][5] 会算出 5， 而它其实应该是 0），所以实现时必须先判断 |i − j| ≤ 1。

反向映射（由 k 求 i, j）。把 k = 2i + j 与 j − i ∈ {−1, 0, 1} 联立，可以得到一个漂亮的闭式：

验证：k = 0 → i = 0, j = 0；k = 2 → i = 1, j = 0； k = 6 → i = 2, j = 2；k = 12 → i = 4, j = 4 ✅ 这条公式在「已知某个一维下标，反求它在矩阵中的位置」这类题里非常好用。

6.6.3 两道手算例题

例 6-10 n = 5 的三对角矩阵按行优先压进 sa[]。 求 a[0][1]、a[2][2]、a[4][3] 的位置， 以及 a[3][0] 的值和数组长度。

例 6-11（反向题） n = 6 的三对角矩阵压缩后得到 sa[0..15]。 问 sa[9] 和 sa[15] 分别对应矩阵中的哪个元素？总长度是多少？

6.6.4 C++ 实现，以及「为什么还能 O(1)」

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   三对角（带状）矩阵的压缩存储（竞赛写法：全局数组 + 映射函数）
   定义：|i - j| <= 1 的元素才可能非零，其余全为 0。
        第 0 行 2 个元素，第 1 .. n-2 行各 3 个，第 n-1 行 2 个，
        总长度 = 2 + 3(n-2) + 2 = 3n - 2。
   （若题目改成「按对角线存进 3×n 的数组」，长度就是 3n，看清题面）

   按行优先压进 sa[0 .. 3n-3]：
       i = 0  : k = j                     （j = 0, 1）
       i >= 1 : k = 3i - 1 + (j - i + 1)  （前 i 行共 2 + 3(i-1) = 3i-1 个元素，
                                            本行内偏移是 j - (i-1)）
   两个式子可以统一写成：k = 2i + j   （在 |i - j| <= 1 的前提下）
   反向公式：i = floor((k+1)/3)，j = k - 2i
   ============================================================ */

const int N = 1005;            /* 最大阶数 */
int sa[3 * N];                 /* 压缩数组，长度 3n - 2（开 3N 个单元够用） */
int n;                         /* 矩阵阶数 */

/* 一维下标公式：k = 2i + j，只要求 |i-j| <= 1 */
inline int idx(int i, int j) { return 2 * i + j; }

/* 带内判断：|i - j| > 1 的位置压缩后根本没存，逻辑上恒为 0 */
inline bool inBand(int i, int j) { return abs(i - j) <= 1; }

/* 反向映射：已知一维下标 k，反求它在矩阵中的位置 (i, j) */
inline void decode(int k, int &i, int &j) { i = (k + 1) / 3; j = k - 2 * i; }

/* 边界检查：竞赛里不抛异常，越界就打印提示 */
bool inRange(int i, int j) {
    if (0 <= i && i < n && 0 <= j && j < n) return true;
    printf("  [warn] (%d, %d) out of range\n", i, j);
    return false;
}

/* 存：带外的位置不存储，只允许写 0（写别的值说明调用者的逻辑已经错了） */
bool setVal(int i, int j, int v) {
    if (!inRange(i, j)) return false;
    if (!inBand(i, j)) {
        if (v != 0) {
            printf("  [warn] a[%d][%d] 在带外，必须是 0\n", i, j);
            return false;
        }
        return true;                     /* 带外的 0 根本不存储 */
    }
    sa[idx(i, j)] = v;
    return true;
}

/* 取：带外一律是 0，这个判断同样是 O(1) */
int getVal(int i, int j) {
    if (!inRange(i, j)) return 0;
    return inBand(i, j) ? sa[idx(i, j)] : 0;
}

/* 打印整个一维数组，方便与手算结果对照 */
void dump() {
    int len = 3 * n - 2;
    for (int k = 0; k < len; ++k) printf("%d%c", sa[k], k + 1 == len ? '\n' : ' ');
}

int main() {
    n = 5;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            if (inBand(i, j)) setVal(i, j, i * 10 + j);   /* a[i][j] = i*10 + j */

    printf("n = 5, storage size = %d (3n - 2 = %d)\n", 3 * n - 2, 3 * n - 2);  /* 13 */
    printf("a[0][1] = %d, k = %d\n", getVal(0, 1), idx(0, 1));   /* 1, 1 */
    printf("a[2][2] = %d, k = %d\n", getVal(2, 2), idx(2, 2));   /* 22, 6 */
    printf("a[4][3] = %d, k = %d\n", getVal(4, 3), idx(4, 3));   /* 43, 11 */
    printf("a[3][0] = %d  (out of band -> 0, not stored)\n", getVal(3, 0));
    printf("sa[] = ");
    dump();                       /* 0 1 10 11 12 21 22 23 32 33 34 43 44 */

    printf("reverse mapping k -> (i, j):\n  ");
    for (int k = 0; k < 3 * n - 2; ++k) {
        int i, j;
        decode(k, i, j);
        printf("k%d=a[%d][%d]  ", k, i, j);
    }
    printf("\n");

    /* 带外写非 0 值会被拦住 */
    printf("setVal(3, 0, 99) -> %s\n", setVal(3, 0, 99) ? "true" : "false");
    return 0;
}

现在回答导读里埋的那个问题：带状矩阵压缩之后，为什么随机存取能力一点没丢？ 原因可以拆成三条：

1. 映射是闭式的。k = 2i + j 只包含一次乘法和一次加法， 与 n 无关，也与「前面存了多少」无关——不需要扫描、不需要累加、不需要查表。 计算机执行它就是一条 lea 指令的事。
2. 带外的元素不需要「找」，逻辑上直接是 0。压缩存储的前提就是 「没存的位置就是 0」，所以 getVal 只要先判断 |i−j| ≤ 1， 不满足就 return 0。这个判断也是 O(1)。
3. 它和原来的二维数组一样是「连续 + 等长」的。压缩后仍然是一个普通的一维数组， sa[k] 本身还是 O(1)。压缩只是把「二维下标」翻译成「一维下标」， 并没有改变底层存储方式。

对比一下稀疏矩阵就清楚了：稀疏矩阵的非零元位置无法用不等式描述， 所以既写不出闭式 k = f(i,j)，也无法用「带外恒为 0」这种简单规则兜底， 只能老老实实把位置存下来再查找。这就是本书 6.7 节要处理的问题。

6.6.5 动画演示：逐行压平与 k 的计算

下面的动画用一个 6 阶三对角矩阵，逐行、逐元素地演示压平过程。 每一帧都会显示 k 是按「前 i 行总数 = 3i − 1」加上「行内偏移 = j − i + 1」 算出来的，请对照公式 k = 2i + j 体会两者为什么等价：

6.7 稀疏矩阵及其压缩存储

6.7.1 什么是稀疏矩阵：公式法为什么在这里失灵

稀疏矩阵（sparse matrix）是指非零元素个数远小于矩阵元素总数的矩阵。 国内教材通常给出一个量化标准：若非零元个数占元素总数的比例不超过 5% （即 tu / (mu × nu) ≤ 0.05），就认为它是稀疏的。

这个「5%」只是个工程上的经验界，不要当成数学定义去较真。真正要抓住的是两点：

• 稀疏是「比例」概念，不是「规模」概念。一个 1000000 × 1000000 的矩阵， 如果每个位置都非零，它依然是稠密矩阵；反过来，一个 10 × 10 的矩阵只要有 2 个非零元， 按比例说它也是稀疏的（只是这时候压缩反而可能不划算）。
• 稀疏矩阵的价值在于「省空间」和「省时间」。工程上（有限元、图计算、推荐系统） 动辄百万阶的稀疏矩阵，用稠密方式存根本开不出内存；而很多算法（比如稀疏矩阵向量乘） 可以做到只与非零元个数成正比，与矩阵阶数几乎无关。

现在回答关键问题：为什么对称矩阵、三角矩阵、三对角矩阵那套「公式法」在稀疏矩阵上不行？

回忆 6.3.3 的三步法。第一步要「确定存哪些元素」，也就是用下标不等式把非零区域描述出来。 对称矩阵能写成 i ≥ j，三对角能写成 |i − j| ≤ 1， 有了这样的不等式，第三步的「数前面有多少个元素」才能变成等差数列求和， 最终得到闭式的 k = f(i,j)。

而一般稀疏矩阵的零元素是东一个西一个、毫无规律的：它可能出现在 (0,3)，也可能出现在 (17,92)，两个非零元之间没有任何位置上的约束关系。 既然写不出不等式，也就数不出「前面有多少个」，公式法直接失灵。

6.7.2 三元组顺序表：最朴素也最常用的方案

三元组顺序表（sequential list of triples）是稀疏矩阵最经典的顺序存储结构。 它的定义短得可以背下来：用一个结构体数组，按行优先顺序 （行号升序；同一行内列号升序）存放全部非零元，每个元素包含 row、col、value 三个域； 另外用 mu、nu、tu 记录行数、列数和非零元个数。

为什么要规定「按行优先有序」？因为有序带来三个直接好处： 可以按行处理（矩阵加法、乘法都要按行归并）、转置时能一次定位、 打印和调试时输出可读。代价是插入或修改一个元素时可能要搬移后续元素—— 和顺序表一样的老毛病。

结构定义与「朴素转置」的实现在下面。请特别留意朴素转置的复杂度是怎么来的—— 它其实就是「按列扫描 + 全表查找」这个笨办法的代价：

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   稀疏矩阵的三元组顺序表 + 朴素转置（竞赛写法：全局数组 + 自由函数）
   三元组：(row, col, value)，按"行优先"排列（行号升序，同行按列号升序）。
   struct Triple 只描述"一个非零元长什么样"，属于结点型，保留；
   矩阵本体就是一个三元组数组 + 三个计数器（mu / nu / tu），不必再套一层 struct。

   朴素转置的思路：
       转置后 T 的第 col 行 = M 的第 col 列。
       于是依次枚举 M 的每一列 col = 0, 1, ..., nu-1，
       为了找到该列的所有非零元，每次都要把整个三元组表扫一遍。
       扫描次数 = nu * tu，所以时间复杂度 O(nu * tu)。
   ============================================================ */

struct Triple {                 /* 结点型：一个非零元 = 行号 + 列号 + 值 */
    int r, c, v;
};

const int MAXT = 100;           /* 数组开够：本题非零元个数 ≤ 100（竞赛里按题目上限开） */

Triple t[MAXT];                 /* 转置结果 T 的三元组表 */
int mu, nu, tu;                 /* T 的行数、列数、非零元个数 */

/* 往 T 里追加一个非零元 —— 三元组表的"插入"就是这一步 */
void pushTriple(int r, int c, int v) {
    t[tu].r = r; t[tu].c = c; t[tu].v = v;
    ++tu;
}

void printMatrix(const char* name) {
    cout << name << " : " << mu << " x " << nu << ", tu = " << tu << "\n";
    for (int p = 0; p < tu; ++p)
        cout << "   [" << p << "] (" << t[p].r << ", " << t[p].c
             << ", " << t[p].v << ")\n";
}

/* 朴素转置：把 src（stu 个非零元、smu 行 snu 列）转置后存进全局的 t[]，O(snu * stu) */
void naiveTranspose(const Triple src[], int stu, int smu, int snu) {
    mu = snu; nu = smu;                 /* 行列互换 */
    tu = stu;
    if (stu == 0) return;
    int q = 0;
    for (int col = 0; col < snu; ++col)               /* 枚举原矩阵的每一列 */
        for (int p = 0; p < stu; ++p)                 /* 每列都重扫整个三元组表 */
            if (src[p].c == col) {
                t[q].r = src[p].c;                    /* 行列互换 */
                t[q].c = src[p].r;
                t[q].v = src[p].v;
                ++q;
            }
}

int main() {
    /* 4 行 5 列，6 个非零元；总元素 20 个，稀疏度 30%（教学用小例子，真实场景常 < 5%） */
    Triple M[MAXT];
    int Mtu = 0;
    M[Mtu++] = {0, 1, 3};
    M[Mtu++] = {0, 4, 7};
    M[Mtu++] = {1, 2, 5};
    M[Mtu++] = {2, 0, 4};
    M[Mtu++] = {2, 3, 6};
    M[Mtu++] = {3, 1, 2};

    /* 为了能同框打印两个矩阵，把 M 的内容也放进全局表里，用 printMatrix 输出 */
    mu = 4; nu = 5; tu = 0;
    for (int p = 0; p < Mtu; ++p) pushTriple(M[p].r, M[p].c, M[p].v);
    printMatrix("M");

    naiveTranspose(M, Mtu, 4, 5);
    printMatrix("T = M^T (naive transpose)");
    return 0;
}

6.7.3 转置：朴素算法为什么慢

矩阵转置（transpose）在数学上很简单：T[j][i] = M[i][j]。 对三元组表来说，把每个三元组的 row 和 col 交换一下就行。 但请注意：光交换还不够。交换之后，数组里的顺序变成了「按原来的列优先」， 不再是「按行优先有序」了，而我们的结构约定要求它必须有序，所以还得重新排。

朴素算法的想法是：不去排序，而是按目标顺序直接产生。 转置后 T 的第 col 行来自 M 的第 col 列， 于是我按 col = 0, 1, …, nu−1 的顺序，依次把 M 中所有列号为 col 的三元组找出来、交换行列、追加到 T 里。 因为「先处理列号小的」，产出的 T 自然就是按行有序的。

问题在于「找出来」这一步：为了找第 col 列的所有非零元， 我们必须把整个三元组表从头到尾扫一遍。一共有 nu 列， 每列都扫 tu 个元素，于是：

这个复杂度有多糟？考虑一个 n × n 的稀疏矩阵，若 tu 与 n 同阶 （比如每行平均只有常数个非零元），那 O(nu × tu) = O(n²)—— 和直接把稠密二维数组转置一遍一样慢，等于白折腾。显然需要改进。

6.7.4 快速转置：num[] 与 cpot[] 的递推

快速转置（fast transpose）的核心洞察只有一句话： 如果能事先知道「转置后第 col 行的第一个元素应该放在哪个位置」， 那么只需要扫一遍 M.data，边扫边放，转置就完成了。

为此需要两个辅助数组，长度都是 nu（转置后的行数）：

num[col]

M 中第 col 列的非零元个数。 因为 M 的列就是 T 的行，所以它同时也是 T 中第 col 行的元素个数。

cpot[col]

T 中第 col 行的第一个元素应该放在 T.data 的哪个下标处（cpot = column position）。

num 很好求：扫一遍 M.data，对每个三元组 ++num[它的 col] 即可。 cpot 则用前缀和递推出来：

这两行的意思很直白：第 col 行的起点 = 第 col−1 行的起点 + 第 col−1 行有多少个元素。 就像排队时分座位：前面那队有多少人，我这队就从下一个空位开始坐。

用图 6-8 里的矩阵 M 实际推一遍（mu=4, nu=5, tu=6）：

• cpot[0] = 0（人为约定：第一行从 0 号位开始）。
• cpot[1] = cpot[0] + num[0] = 0 + 1 = 1。
• cpot[2] = cpot[1] + num[1] = 1 + 2 = 3。
• cpot[3] = cpot[2] + num[2] = 3 + 1 = 4。
• cpot[4] = cpot[3] + num[3] = 4 + 1 = 5。
• 校验：cpot[4] + num[4] = 5 + 1 = 6 = tu ✅ 刚好把 6 个位置分配完， 这也是一条常用的自检规则。

快速转置的完整实现如下，三个步骤在注释里标了序号：

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   稀疏矩阵的快速转置 O(nu + tu)（竞赛写法：全局数组 + 自由函数）
   核心：先算出"转置后每一行从哪里开始放"，就能一次扫描直接定位。

   num[col]  = 原矩阵第 col 列的非零元个数（也就是 T 中第 col 行的元素个数）
   cpot[col] = T 中第 col 行第一个元素在 t[] 中的下标（cpot = column position）
       递推：cpot[0] = 0
             cpot[col] = cpot[col-1] + num[col-1]     (col = 1 .. nu-1)

   然后按原三元组表的顺序扫一遍：
       对每个三元组 (r, c, v)，它应该放在 t[cpot[c]]，
       放完后 ++cpot[c]，因为这一行的下一个元素要往后挪一格。
   ============================================================ */

struct Triple { int r, c, v; };   /* 结点型：一个非零元 */

const int MAXT = 100;             /* 数组开够：本题非零元个数 ≤ 100 */
const int MAXC = 100;             /* 列数上限：num / cpot 的长度依据 */

Triple t[MAXT];                   /* 转置结果 T 的三元组表 */
int mu, nu, tu;                   /* T 的行数、列数、非零元个数 */
int num[MAXC], cpot[MAXC];        /* ① 每列非零元个数  ② 每行起始下标 */

void pushTriple(int r, int c, int v) {   /* 三元组表的插入操作 */
    t[tu].r = r; t[tu].c = c; t[tu].v = v;
    ++tu;
}

void printMatrix(const char* name) {
    cout << name << " : " << mu << " x " << nu << ", tu = " << tu << "\n";
    for (int p = 0; p < tu; ++p)
        cout << "   [" << p << "] (" << t[p].r << ", " << t[p].c << ", " << t[p].v << ")\n";
}

/* 就地求 num[] / cpot[]：把 src 的列统计出来，供按行定位使用 */
void buildNumCpot(const Triple src[], int stu, int snu) {
    for (int col = 0; col < snu; ++col) num[col] = 0;
    for (int p = 0; p < stu; ++p) ++num[src[p].c];          /* ① 统计每列的非零元个数 */

    if (snu > 0) cpot[0] = 0;
    for (int col = 1; col < snu; ++col)                     /* ② 递推每行的起始位置 */
        cpot[col] = cpot[col - 1] + num[col - 1];
}

void fastTranspose(const Triple src[], int stu, int smu, int snu, bool verbose = false) {
    mu = snu; nu = smu;                                     /* 行列互换 */
    tu = stu;
    if (stu == 0) return;

    buildNumCpot(src, stu, snu);
    if (verbose) {
        cout << "col :"; for (int c = 0; c < snu; ++c) cout << "  " << c; cout << "\n";
        cout << "num :"; for (int c = 0; c < snu; ++c) cout << "  " << num[c]; cout << "\n";
        cout << "cpot:"; for (int c = 0; c < snu; ++c) cout << "  " << cpot[c]; cout << "\n";
    }

    for (int p = 0; p < stu; ++p) {                         /* ③ 一次扫描，直接定位 */
        int col = src[p].c;
        int q = cpot[col];
        t[q].r = src[p].c;
        t[q].c = src[p].r;
        t[q].v = src[p].v;
        ++cpot[col];                                        /* 该行的下一个空位 */
    }
}

int main() {
    Triple M[MAXT];
    int Mtu = 0;
    M[Mtu++] = {0, 1, 3};
    M[Mtu++] = {0, 4, 7};
    M[Mtu++] = {1, 2, 5};
    M[Mtu++] = {2, 0, 4};
    M[Mtu++] = {2, 3, 6};
    M[Mtu++] = {3, 1, 2};

    mu = 4; nu = 5; tu = 0;
    for (int p = 0; p < Mtu; ++p) pushTriple(M[p].r, M[p].c, M[p].v);
    printMatrix("M");

    fastTranspose(M, Mtu, 4, 5, true);
    printMatrix("T = M^T (fast transpose)");

    /* 快速转置的结果与"按行扫描"的期望完全一致：
       T 的第 0 行来自 M 的第 0 列 ... 依次类推 */
    return 0;
}

6.7.5 行逻辑链接的顺序表：给三元组表加一个 rpos[]

三元组表虽然按行有序，但「第 i 行从哪儿开始」这件事它没说。 想知道第 3 行有哪些元素，你还是得从头扫或者二分查找，最坏 O(tu)。 解决办法和快速转置里用的一模一样——再加一个前缀和数组：

建立了 rpos 之后，第 i 行的所有非零元就是 data[rpos[i] .. rpos[i+1]−1]，访问代价降为 O(该行非零元个数)。这种带行起始位置的三元组表，教材里叫 行逻辑链接的顺序表（row-logical linked list）—— 名字里的「链接」指的就是这个 rpos 起到的「逻辑链接」作用， 它并不是真的指针。

有了它，稀疏矩阵加法就能像归并两个有序数组那样按行做了： 对每一行 i，用两个指针分别扫 A 和 B 在该行的元素， 列号小的先输出、列号相等就相加。复杂度 O(tuA + tuB)。

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   行逻辑链接的顺序表（Row-Logical Linked List）
   在三元组顺序表上增加一个 rpos[] 数组：
       rpos[i] = 第 i 行第一个非零元在 t[] 中的下标，rpos[mu] = tu
   递推：先数出每行的非零元个数，再做前缀和。
   有了 rpos[]，遍历第 i 行只需 O(该行非零元个数)，不用再扫全表。
   有了它，稀疏矩阵加法就能按行归并：O(tuA + tuB)。

   竞赛写法：A、B、C 三个矩阵各用"一组全局数组 + 三个计数器"表示，
   不再有 SparseMatrix 这层 struct 壳；struct Triple 描述非零元，保留。
   ============================================================ */

struct Triple { int r, c, v; };   /* 结点型：一个非零元 = 行号 + 列号 + 值 */

const int MAXT = 100;             /* 数组开够：每个矩阵非零元个数 ≤ 100 */
const int MAXR = 100;             /* 行数上限：rpos 的长度依据 */

/* ---------- A ---------- */
Triple ta[MAXT];
int muA, nuA, tuA;
int rposA[MAXR];
/* ---------- B ---------- */
Triple tb[MAXT];
int muB, nuB, tuB;
int rposB[MAXR];
/* ---------- C = A + B ---------- */
Triple tc[MAXT];
int muC, nuC, tuC;
int rposC[MAXR];

void pushA(int r, int c, int v) { ta[tuA].r = r; ta[tuA].c = c; ta[tuA].v = v; ++tuA; }
void pushB(int r, int c, int v) { tb[tuB].r = r; tb[tuB].c = c; tb[tuB].v = v; ++tuB; }
void pushC(int r, int c, int v) { tc[tuC].r = r; tc[tuC].c = c; tc[tuC].v = v; ++tuC; }

/* 建 rpos[]：两遍扫描，O(mu + tu) */
void buildRposA() {
    for (int i = 0; i <= muA; ++i) rposA[i] = 0;
    for (int p = 0; p < tuA; ++p) ++rposA[ta[p].r + 1];       /* ① 数每行的个数 */
    for (int i = 1; i <= muA; ++i) rposA[i] += rposA[i - 1];  /* ② 前缀和 */
}
void buildRposB() {
    for (int i = 0; i <= muB; ++i) rposB[i] = 0;
    for (int p = 0; p < tuB; ++p) ++rposB[tb[p].r + 1];
    for (int i = 1; i <= muB; ++i) rposB[i] += rposB[i - 1];
}
void buildRposC() {
    for (int i = 0; i <= muC; ++i) rposC[i] = 0;
    for (int p = 0; p < tuC; ++p) ++rposC[tc[p].r + 1];
    for (int i = 1; i <= muC; ++i) rposC[i] += rposC[i - 1];
}

void printA(const char* name) {
    cout << name << " : " << muA << " x " << nuA << ", tu = " << tuA << "\n";
    for (int p = 0; p < tuA; ++p)
        cout << "   [" << p << "] (" << ta[p].r << ", " << ta[p].c << ", " << ta[p].v << ")\n";
    cout << "   rpos =";
    for (int i = 0; i <= muA; ++i) cout << " " << rposA[i];
    cout << "\n";
}
void printB(const char* name) {
    cout << name << " : " << muB << " x " << nuB << ", tu = " << tuB << "\n";
    for (int p = 0; p < tuB; ++p)
        cout << "   [" << p << "] (" << tb[p].r << ", " << tb[p].c << ", " << tb[p].v << ")\n";
    cout << "   rpos =";
    for (int i = 0; i <= muB; ++i) cout << " " << rposB[i];
    cout << "\n";
}
void printC(const char* name) {
    cout << name << " : " << muC << " x " << nuC << ", tu = " << tuC << "\n";
    for (int p = 0; p < tuC; ++p)
        cout << "   [" << p << "] (" << tc[p].r << ", " << tc[p].c << ", " << tc[p].v << ")\n";
    cout << "   rpos =";
    for (int i = 0; i <= muC; ++i) cout << " " << rposC[i];
    cout << "\n";
}

/* 按行访问 A：打印第 i 行的所有非零元，复杂度 O(该行非零元个数) */
void printRowA(int i) {
    cout << "row " << i << " :";
    for (int p = rposA[i]; p < rposA[i + 1]; ++p)
        cout << "  col" << ta[p].c << "=" << ta[p].v;
    cout << "\n";
}

/* 稀疏矩阵加法 C = A + B（同型），按行归并，O(tuA + tuB) */
void addMatrix() {
    if (muA != muB || nuA != nuB) { printf("shape mismatch\n"); return; }
    muC = muA; nuC = nuA; tuC = 0;
    for (int i = 0; i < muA; ++i) {
        int pa = rposA[i], ea = rposA[i + 1];     /* A 的第 i 行区间 */
        int pb = rposB[i], eb = rposB[i + 1];     /* B 的第 i 行区间 */
        while (pa < ea && pb < eb) {
            int ca = ta[pa].c, cb = tb[pb].c;
            if (ca < cb) { pushC(i, ca, ta[pa].v); ++pa; }
            else if (ca > cb) { pushC(i, cb, tb[pb].v); ++pb; }
            else {
                int s = ta[pa].v + tb[pb].v;
                if (s != 0) pushC(i, ca, s);   /* 相加为 0 就不必存了 —— 稀疏矩阵加法的关键细节 */
                ++pa; ++pb;
            }
        }
        for (; pa < ea; ++pa) pushC(i, ta[pa].c, ta[pa].v);
        for (; pb < eb; ++pb) pushC(i, tb[pb].c, tb[pb].v);
    }
}

int main() {
    muA = 4; nuA = 5; tuA = 0;
    pushA(0, 1, 3); pushA(0, 4, 7); pushA(1, 2, 5);
    pushA(2, 0, 4); pushA(2, 3, 6); pushA(3, 1, 2);

    muB = 4; nuB = 5; tuB = 0;
    pushB(0, 4, -7); pushB(1, 2, 5); pushB(2, 0, -4);
    pushB(2, 3, 6);  pushB(3, 2, 9);

    buildRposA(); buildRposB();
    printA("A"); printB("B");

    printRowA(0);
    printRowA(2);

    addMatrix();
    buildRposC();
    printC("C = A + B");

    /* 注意看：C 的第 0 行只剩 (0,1,3)，
       因为 A(0,4,7) + B(0,4,-7) = 0，这个位置直接"消失"了。 */
    return 0;
}

运行结果里最值得看的一处是：C 的第 0 行只剩下 (0,1,3) 一项， 因为 A(0,4,7) + B(0,4,−7) = 0，这个位置在结果里直接消失了。 代码里那一句 if (s != 0) C.push(...) 就是稀疏矩阵运算的典型风格： 结果为零的位置绝不存储。

6.7.6 十字链表：能按行也能按列，还便于增删

三元组表有两个绕不开的弱点：一是插入删除要搬移元素，二是按列访问很不方便。 而矩阵运算（尤其是乘法、加法这类会「边算边生成新非零元」的场景）恰恰需要频繁插入， 也常常需要按列扫描。于是有了链式方案——十字链表（orthogonal list）。

它把矩阵的每个非零元做成一个结点，结点有五个域： row（行号）、col（列号）、value（值）、 down（指向同列下一个非零元）、right（指向同行下一个非零元）。 再加上两个表头指针数组：

• rhead[mu]：rhead[i] 指向第 i 行的第一个非零元， 顺着 right 就能走完整行；
• chead[nu]：chead[j] 指向第 j 列的第一个非零元， 顺着 down 就能走完整列。

两条链表在同一个结点上「十字交叉」，所以叫十字链表。 它的最大好处是：插入或删除一个非零元时，只需要改动常数个指针，不用搬移任何数据 （这也是所有链式结构的共同优点），而且「按行」和「按列」两条访问路径都是现成的。

实现十字链表的关键技巧是「用指针的指针找插入点」： 要往一条有序链表里插入，不必区分「插在头部」和「插在中间」两种情况， 只要用一个二级指针 OLNode** p 一路往后走，最后 *p = node 就完事了。 下面的代码同时演示了「按行列有序插入」和「同一位置重复插入时累加」两个细节：

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   十字链表（Orthogonal List）—— 稀疏矩阵的链式存储（竞赛写法）
   每个非零元是一个含 5 个域的结点：
       row, col, value, down, right
   right 把同一行的非零元串成链表（按列号升序）；
   down  把同一列的非零元串成链表（按行号升序）。
   另外用两个指针数组作为表头：rhead[mu]（行）、chead[nu]（列）。

   优点：插入 / 删除非零元不用像顺序表那样成片搬移元素，
         而且「按行」和「按列」两条访问路径都是 O(1) 起步。
   代价：每个结点多两个指针域；按下标 (i,j) 随机取值退化成沿链查找。
   ============================================================ */

const int N = 1005;                 /* 行数、列数上限 */
struct OLNode {                     /* 结点：5 个域 */
    int r, c, v;                    /* 行号、列号、值 */
    OLNode *down, *right;           /* down = 同列下一个，right = 同行下一个 */
};

OLNode *rhead[N], *chead[N];        /* 行 / 列表头指针数组 */
int mu, nu, tu;                     /* 行数、列数、非零元个数 */

/* 初始化成 mu 行 nu 列的空矩阵：表头全部置空 */
void initMatrix(int m, int n) {
    mu = m; nu = n; tu = 0;
    for (int i = 0; i < mu; ++i) rhead[i] = nullptr;
    for (int j = 0; j < nu; ++j) chead[j] = nullptr;
}

/* 新建结点；竞赛里干脆不写 delete，程序结束内存就还给操作系统了 */
OLNode* newNode(int r, int c, int v) {
    OLNode* p = new OLNode;
    p->r = r; p->c = c; p->v = v;
    p->down = p->right = nullptr;
    return p;
}

/* 插入（或累加）一个非零元，保持行、列两条链表都有序 */
void insertNode(int r, int c, int v) {
    if (v == 0) return;                        /* 0 不是非零元，不入表 */

    /* ① 先在第 r 行的行链表里找插入点（按列号升序）
          「用指针的指针找插入点」：不必区分插在头部还是中间，
          一路往后走，最后 *rp = node 一句就搞定 */
    OLNode** rp = &rhead[r];
    while (*rp && (*rp)->c < c) rp = &(*rp)->right;
    if (*rp && (*rp)->c == c) {                /* 该位置已经有结点了 → 只累加值 */
        (*rp)->v += v;                         /* 列链表里挂的是同一个结点，不用再动 */
        return;
    }
    OLNode* node = newNode(r, c, v);
    node->right = *rp;                         /* 先接后面，再改前面，顺序不能反 */
    *rp = node;

    /* ② 再挂到第 c 列的列链表（按行号升序） */
    OLNode** cp = &chead[c];
    while (*cp && (*cp)->r < r) cp = &(*cp)->down;
    node->down = *cp;
    *cp = node;

    ++tu;
}

void printByRow() {
    printf("by row:\n");
    for (int i = 0; i < mu; ++i) {
        printf("  row %d :", i);
        for (OLNode* p = rhead[i]; p; p = p->right)
            printf("  (%d,%d,%d)", p->r, p->c, p->v);
        printf("\n");
    }
}
void printByCol() {
    printf("by col:\n");
    for (int j = 0; j < nu; ++j) {
        printf("  col %d :", j);
        for (OLNode* p = chead[j]; p; p = p->down)
            printf("  (%d,%d,%d)", p->r, p->c, p->v);
        printf("\n");
    }
}
int nonzeros() { return tu; }

int main() {
    initMatrix(4, 5);
    insertNode(0, 1, 3);
    insertNode(0, 4, 7);
    insertNode(1, 2, 5);
    insertNode(2, 0, 4);
    insertNode(2, 3, 6);
    insertNode(3, 1, 2);
    insertNode(0, 1, 10);          /* 同一位置的再次插入 → 值累加为 13 */
    insertNode(1, 1, 1);           /* 插在已存在的列链表中间，链式结构无需搬移 */

    printf("nonzeros = %d\n", nonzeros());
    printByRow();
    printByCol();
    return 0;
}

6.7.7 工程视角：CSR / CSC 格式

前面讲的三种结构都是教材语言。真到了工程里（SciPy、PETSc、cuSPARSE、PyTorch 的稀疏张量）， 用得最多的是 CSR（Compressed Sparse Row，压缩稀疏行）格式。 它其实是「三元组表 + rpos」的极致版本——把行号也省掉了（三个数组与原矩阵的对应关系见图 6-11）：

注意 rowPtr 和行逻辑链接里的 rpos 是同一个东西， 只是多存了一个哨兵 rowPtr[m] = nnz，这样每行的区间都能写成半开区间 [rowPtr[i], rowPtr[i+1])，代码里少一个 if。

为什么说 CSR 是 cache 友好的？看稀疏矩阵向量乘 y = Ax 的循环：

内层循环在 values[] 和 colIdx[] 上是严格顺序前进的， 这正是 CPU 预取器和缓存最喜欢的访问模式：读入 values[p] 时， 同一缓存行里的后面十几个元素也一起进了缓存，接下来的迭代几乎是零开销。 相比之下，十字链表每访问一个元素都要解引用一个指针，地址完全随机， 每次都可能是一次缓存未命中（cache miss）；而三元组表虽然也顺序， 但每个元素多存了一个行号，且没有「按行定位」的能力。

代价是 CSR 的「结构是静态的」：rowPtr 靠前缀和一次性算好， 中途插入一个非零元会导致后面所有位置整体后移。所以 CSR 适合「结构定下来之后反复做运算」 （迭代求解、稀疏矩阵乘法、图算法），不适合「边算边改结构」。 列优先的版本叫 CSC（Compressed Sparse Column），把行、列的角色对调， 适用于按列访问更频繁的场合（比如 MATLAB 历史上就用 CSC）。

#include <iostream>
#include <vector>
#include <iomanip>
using namespace std;

/* ============================================================
   CSR（Compressed Sparse Row，压缩稀疏行）
   现代工程里最常用的稀疏矩阵格式，由三个一维数组构成：

       values[] : 按行优先顺序排列的全部非零元的值
       colIdx[] : 与 values[] 一一对应的列号
       rowPtr[] : 长度为 m+1；第 i 行的非零元在 values[] 中的区间是
                  [rowPtr[i], rowPtr[i+1])

   为什么 cache 友好：
       第 i 行的数据在 values[]/colIdx[] 里是"连续的一整段"，
       遍历一行时访存完全顺序，不像十字链表那样满内存跳指针；
       而且除了这三个数组外没有任何指针开销。
   ============================================================ */

struct Triple { int r, c, v; };

struct CSR {
    int m, n;
    vector<int> values, colIdx, rowPtr;
};

CSR buildCSR(int m, int n, const vector<Triple>& trips) {
    CSR A; A.m = m; A.n = n;
    A.rowPtr.assign(m + 1, 0);
    for (size_t p = 0; p < trips.size(); ++p) ++A.rowPtr[trips[p].r + 1];   /* ① 数每行个数 */
    for (int i = 1; i <= m; ++i) A.rowPtr[i] += A.rowPtr[i - 1];           /* ② 前缀和 */

    A.values.assign(trips.size(), 0);
    A.colIdx.assign(trips.size(), 0);
    vector<int> pos(A.rowPtr.begin(), A.rowPtr.end() - 1);                 /* ③ 每行的写入位置 */
    for (size_t p = 0; p < trips.size(); ++p) {
        int q = pos[trips[p].r]++;
        A.values[q] = trips[p].v;
        A.colIdx[q] = trips[p].c;
    }
    return A;
}

void printCSR(const CSR& A, const char* name) {
    cout << name << " (" << A.m << " x " << A.n << ", nnz = "
         << A.values.size() << ")\n";
    cout << "  rowPtr ="; for (size_t i = 0; i < A.rowPtr.size(); ++i) cout << " " << A.rowPtr[i]; cout << "\n";
    cout << "  colIdx ="; for (size_t i = 0; i < A.colIdx.size(); ++i) cout << " " << A.colIdx[i]; cout << "\n";
    cout << "  values ="; for (size_t i = 0; i < A.values.size(); ++i) cout << " " << A.values[i]; cout << "\n";
}

/* 稀疏矩阵 × 稠密向量 y = A x，复杂度 O(nnz)，是迭代法解线性方程组的核心算子 */
vector<double> spmv(const CSR& A, const vector<double>& x) {
    vector<double> y(A.m, 0.0);
    for (int i = 0; i < A.m; ++i)
        for (int p = A.rowPtr[i]; p < A.rowPtr[i + 1]; ++p)
            y[i] += A.values[p] * x[A.colIdx[p]];
    return y;
}

int main() {
    vector<Triple> trips = {
        {0, 1, 3}, {0, 4, 7}, {1, 2, 5},
        {2, 0, 4}, {2, 3, 6}, {3, 1, 2}
    };
    CSR A = buildCSR(4, 5, trips);
    printCSR(A, "A");

    vector<double> x = {1, 2, 3, 4, 5};
    vector<double> y = spmv(A, x);
    cout << "y = A x =";
    for (size_t i = 0; i < y.size(); ++i) cout << " " << y[i];
    cout << "\n";     /* 41 15 28 4 —— 手算：3*2+7*5=41，5*3=15，4*1+6*4=28，2*2=4 */
    return 0;
}

6.7.8 稀疏矩阵的加法与乘法

加法在 6.7.5 已经见过：按行归并，O(tuA + tuB)， 关键是「结果为零的位置不存」。这里再强调一个实战中的坑： 稀疏 + 稀疏未必还是稀疏。如果两个矩阵的非零位置几乎不重叠， 结果的非零元个数可能是两者之和（这叫 fill-in，填充）， 极端情况下结果会变成稠密矩阵。所以工程库在处理前通常会先估算结果的 nnz 上界。

乘法稍微绕一点。设 C = A × B，A 是 m × n， B 是 n × p。如果直接照抄三重循环的稠密写法， 复杂度是 O(m × n × p)，完全没利用稀疏性。 正确姿势是「以 A 的非零元为驱动」：

1. 开一个长度为 p 的累加器数组 ctemp[]，初值全 0；
2. 对 A 的第 i 行，遍历它的每个非零元 a[i][k]；
3. 把 a[i][k] 乘到 B 第 k 行的每个非零元 b[k][j] 上， 累加进 ctemp[j]；
4. 第 i 行处理完，把 ctemp 中的非零项收进 C 的第 i 行，然后把 ctemp 清零。

复杂度是 O(m × p + Σa[i][k]≠0 (B 第 k 行的非零元个数))—— 只有当 B 的行很稠密时才会退化。这一段代码里 rpos[] 再次出场， 用来 O(1) 定位 B 的每一行：

#include <bits/stdc++.h>
using namespace std;

/* ============================================================
   稀疏矩阵乘法 C = A x B（行逻辑链接版，竞赛写法：全局数组 + 自由函数）
   设 A 是 m x n，B 是 n x p。
   思路：对 A 的第 i 行的每个非零元 a[i][k]，
         把它乘到 B 的第 k 行的每个非零元 b[k][j] 上，
         结果累加到临时行 ctemp[j] 里；第 i 行处理完，
         把 ctemp 中的非零项收集成 C 的第 i 行。
   复杂度：O(m*p + sum_{a[i][k] != 0} (B 第 k 行的非零元个数))

   这里 A、B、C 各用"三元组全局数组 + rpos[] + 三个计数器"表示；
   struct Triple 只描述非零元，属于结点型，保留。
   ============================================================ */

struct Triple { int r, c, v; };   /* 结点型：一个非零元 */

const int MAXT = 100;             /* 数组开够：每个矩阵非零元个数 ≤ 100 */
const int MAXR = 100;             /* 行数上限：rpos 的长度依据 */
const int MAXC = 100;             /* 列数上限：ctemp 的长度依据 */

/* ---------- A ---------- */
Triple ta[MAXT];
int muA, nuA, tuA, rposA[MAXR];
/* ---------- B ---------- */
Triple tb[MAXT];
int muB, nuB, tuB, rposB[MAXR];
/* ---------- C = A x B ---------- */
Triple tc[MAXT];
int muC, nuC, tuC, rposC[MAXR];

int ctemp[MAXC];                  /* 累加 C 当前这一行的结果 */

void pushA(int r, int c, int v) { ta[tuA].r = r; ta[tuA].c = c; ta[tuA].v = v; ++tuA; }
void pushB(int r, int c, int v) { tb[tuB].r = r; tb[tuB].c = c; tb[tuB].v = v; ++tuB; }
void pushC(int r, int c, int v) { tc[tuC].r = r; tc[tuC].c = c; tc[tuC].v = v; ++tuC; }

void buildRposA() {
    for (int i = 0; i <= muA; ++i) rposA[i] = 0;
    for (int p = 0; p < tuA; ++p) ++rposA[ta[p].r + 1];
    for (int i = 1; i <= muA; ++i) rposA[i] += rposA[i - 1];
}
void buildRposB() {
    for (int i = 0; i <= muB; ++i) rposB[i] = 0;
    for (int p = 0; p < tuB; ++p) ++rposB[tb[p].r + 1];
    for (int i = 1; i <= muB; ++i) rposB[i] += rposB[i - 1];
}
void buildRposC() {
    for (int i = 0; i <= muC; ++i) rposC[i] = 0;
    for (int p = 0; p < tuC; ++p) ++rposC[tc[p].r + 1];
    for (int i = 1; i <= muC; ++i) rposC[i] += rposC[i - 1];
}

void printA(const char* name) {
    cout << name << " : " << muA << " x " << nuA << ", tu = " << tuA << "\n";
    for (int p = 0; p < tuA; ++p)
        cout << "   (" << ta[p].r << ", " << ta[p].c << ", " << ta[p].v << ")\n";
}
void printB(const char* name) {
    cout << name << " : " << muB << " x " << nuB << ", tu = " << tuB << "\n";
    for (int p = 0; p < tuB; ++p)
        cout << "   (" << tb[p].r << ", " << tb[p].c << ", " << tb[p].v << ")\n";
}
void printC(const char* name) {
    cout << name << " : " << muC << " x " << nuC << ", tu = " << tuC << "\n";
    for (int p = 0; p < tuC; ++p)
        cout << "   (" << tc[p].r << ", " << tc[p].c << ", " << tc[p].v << ")\n";
}

void multiply() {
    if (nuA != muB) { printf("shape mismatch: A.nu != B.mu\n"); return; }
    muC = muA; nuC = nuB; tuC = 0;
    if (tuA == 0 || tuB == 0) return;

    for (int i = 0; i < muA; ++i) {
        for (int j = 0; j < nuB; ++j) ctemp[j] = 0;    /* 清空累加器 */
        for (int pa = rposA[i]; pa < rposA[i + 1]; ++pa) {
            int k = ta[pa].c;                         /* A 的列 = B 的行 */
            for (int pb = rposB[k]; pb < rposB[k + 1]; ++pb) {
                int j = tb[pb].c;
                ctemp[j] += ta[pa].v * tb[pb].v;
            }
        }
        for (int j = 0; j < nuB; ++j)
            if (ctemp[j] != 0) pushC(i, j, ctemp[j]);  /* 只保留非零元 */
    }
    buildRposC();
}

int main() {
    muA = 3; nuA = 4; tuA = 0;
    pushA(0, 0, 1); pushA(0, 2, 2); pushA(1, 1, 3);
    pushA(2, 0, 4); pushA(2, 3, 5);

    muB = 4; nuB = 2; tuB = 0;
    pushB(0, 0, 1); pushB(0, 1, 2); pushB(1, 1, 1);
    pushB(2, 0, 3); pushB(3, 0, 1); pushB(3, 1, 1);

    buildRposA(); buildRposB();
    printA("A"); printB("B");

    multiply();
    printC("C = A x B");
    /* 手算对照：
       C[0][0] = 1*1 + 2*3 = 7      C[0][1] = 1*2 + 2*0 = 2
       C[1][0] = 3*0 = 0（不存）    C[1][1] = 3*1 = 3
       C[2][0] = 4*1 + 5*1 = 9      C[2][1] = 4*2 + 5*1 = 13  */
    return 0;
}

6.7.9 动画演示：稀疏矩阵到三元组表的转换

最后用一个动画把「扫描稀疏矩阵、生成三元组表」的过程走一遍。 请留意两点：扫描顺序是行优先（保证了结果有序）， 以及零元素被直接跳过（这就是省空间的来源）：

6.8 四种特殊矩阵横向对比

把本章讲过的方案放进一张表里对照。这张表建议直接背下来—— 考试里的计算题、选择题几乎都能在里面找到答案：

6.9 工程视角：矩阵压缩在数值计算与推荐系统里怎么用

前面八节把矩阵压缩当成「考试计算题」讲完了。但矩阵和稀疏矩阵在真实世界里是另一种身份： 它们是科学计算、机器学习、图算法的基础设施。NumPy / SciPy 的底层是 BLAS 与 LAPACK， GPU 上的 cuSPARSE 与 cuBLAS，PageRank 的一次迭代、有限元的一次求解、推荐系统的一次召回， 全都落在本章讲的这些存储结构上。这一节换个视角，对六个工程落点逐条回答三个问题： 用什么结构 → 为什么必须用它 → 代价和陷阱是什么。

6.9.1 行优先 vs 列优先：一次遍历凭什么差几倍

6.2 节讲了两种压平方式只是地址公式不同。到了工程里，这个「不同」直接换成钱和时间。 先看两个数量级（选 double，每个 8 字节）：

现在拿一个 1024 × 1024 的 double 矩阵做实验。 它占 1024 × 1024 × 8 = 8388608 字节，也就是 8 MB， 远超典型 256 KB ~ 1 MB 的 L2 缓存，整个矩阵只能待在主存里。两种遍历方式写的代码几乎一样：

按行遍历时，8 个 double 挤在同一条缓存行里，平均每 8 次访问才发生一次主存访问； 按列遍历时，a[0][j] 与 a[1][j] 在内存里正好相隔 8 KB， 每一次访问都落在不同的缓存行上，于是 100 万次访问里有 100 万次是主存访问。 把 n = 1024 代进上面两个数字：1024 × 1024 ÷ 8 = 131072 次缓存未命中， 每次 60 ~ 100 ns，光等内存就是 8 ~ 13 ms；而按行遍历只有约 16 万次「有效访问」、 其中绝大多数命中 L1，总耗时不到 1 个周期 × 100 万次的量级。 实测结论很稳定：同一个矩阵、同一个求和结果，按错误顺序遍历会慢 3 ~ 5 倍。

这跟第 01 讲讲的时间复杂度并不矛盾：两种写法的抽象操作次数完全一样， 变的只是每次操作的常数因子，而常数因子里藏着 18 倍的缓存代价。 大 O 描述增长趋势，缓存决定实际快慢，工程优化往往就发生在后者。

回到语言层面，这个坑之所以常见，是因为「哪一维连续」并不统一： C / C++ 的行优先把一行放在一起，而 Fortran、MATLAB、以及 NumPy 的默认布局是列优先。 同一个 A[3][4]，C 里 A[1][2] 落在 (1×4 + 2)×L， Fortran 里落在 (2×3 + 1)×L——不是同一个位置。

NumPy 把这件事显式暴露出来：np.array(a, order='C') 是行优先（C order）， np.asfortranarray 与 order='F' 是列优先（Fortran order，F order）。 同一个矩阵两种布局，内存占用一模一样，但「按行求和」与「按列求和」谁快， 结果会整个反过来。检查一个数组到底是什么布局，靠 .flags['C_CONTIGUOUS']， 或者干脆看 .strides——步长是理解布局最直接的方式， 它比记「行优先还是列优先」更不容易搞错。

还有一条更真实的工程经验：调 BLAS 时，「矩阵该不该转置」是一个性能问题，而不是风格问题。 BLAS 为两种语言习惯提供了两套入口（dgemm 与 CBLAS 的 CblasRowMajor / CblasColMajor），LAPACK 则只有列优先版本。于是：

• 按底层布局调库：需要行优先入口时用 CBLAS 的 RowMajor 包装， 代价是库内部可能再做一次隐式拷贝。
• 转置成库要的布局：拷贝成本是 2mn 次读写。 一个 5000 × 5000 的矩阵要搬 2.5 × 10⁷ 个 double， 而同样规模的一次矩阵乘法约 2n³ = 2.5 × 1011 次浮点运算—— 单次调用时转置成本只有万分之一，可以忽略。
• 但会被调用次数放大：迭代求解器里 SpMV 要跑几万次， 每次都先转置的话搬运量就与乘法本身的运算量同量级。 这就是「一次性转换格式、之后反复用」的由来。

6.9.2 对称矩阵的工程意义：一半内存，以及一半的麻烦

6.4 节说对称矩阵可以只存下三角，长度从 n² 降到 n(n+1)/2。 在考试里这是省一半空间；在工程里，这条性质经常是「算得动」和「算不动」的分界。 三个典型来源：

• 有限元 / 结构力学的刚度矩阵 K：单元之间的相互作用是对称的 （节点 i 对 j 的作用与 j 对 i 的作用成对出现），所以 K 天生对称。 一个三维实体网格动辄十万到百万个自由度，K 是这里最大的单个数据结构。
• 协方差矩阵：Cov(Xi, Xj) = Cov(Xj, Xi)， 恒对称；PCA 要对它做特征分解。
• 距离矩阵 / 相似度矩阵：d(i,j) = d(j,i) 恒成立， 聚类、最近邻、图嵌入里都要把它一次性算出来。

/* ============================================================
   csr_spmv.cpp —— COO → CSR 转换 + 稀疏矩阵乘向量 SpMV

   输入：5 x 5 稀疏矩阵的三元组（行、列、值）
   输出：CSR 的三个数组 values / colIdx / rowPtr，以及 y = A * x

   编译： g++ -std=c++17 -O2 -o csr_spmv csr_spmv.cpp
   ============================================================ */
#include <cstdio>

const int MAXN   = 100;    /* 行数上界 */
const int MAXNNZ = 1000;   /* 非零元个数上界 */

/* COO 三元组：行、列、值 */
struct Triple { int row, col, val; };

int values[MAXNNZ];      /* 非零元的值 */
int colIdx[MAXNNZ];      /* 与 values 一一对应的列号 */
int rowPtr[MAXN + 1];    /* 第 i 行的区间是 [rowPtr[i], rowPtr[i+1]) */

/* ---------- COO -> CSR：计数 + 前缀和 + 按行放置 ---------- */
void cooToCsr(const Triple t[], int nnz, int rows) {
    for (int i = 0; i <= rows; ++i) rowPtr[i] = 0;

    /* 1) 统计每一行有几个非零元：先记在 rowPtr[row + 1] 上 */
    for (int k = 0; k < nnz; ++k) rowPtr[t[k].row + 1]++;

    /* 2) 前缀和：rowPtr[i + 1] = 0..i 行的非零元总数，即第 i+1 行的起点 */
    for (int i = 0; i < rows; ++i) rowPtr[i + 1] += rowPtr[i];

    /* 3) 再扫一遍三元组，把元素放进它该在的位置。
          rowPtr[row] 先当写入游标，写一个就往右挪一格，
          全部写完时它恰好等于第 row+1 行的起点 */
    for (int k = 0; k < nnz; ++k) {
        int p = rowPtr[t[k].row]++;
        values[p] = t[k].val;
        colIdx[p] = t[k].col;
    }

    /* 4) 把游标整体右移一位还原成行起点，并补上哨兵 rowPtr[rows] = nnz */
    for (int i = rows; i > 0; --i) rowPtr[i] = rowPtr[i - 1];
    rowPtr[0] = 0;
}

/* ---------- SpMV：y = A * x，按行扫描 CSR ---------- */
void spmvCsr(int rows, const int x[], int y[]) {
    for (int i = 0; i < rows; ++i) {
        int sum = 0;
        /* 第 i 行的非零元就是 values[rowPtr[i] .. rowPtr[i+1] - 1] 这一段，
           连续访问，不需要任何查找 */
        for (int p = rowPtr[i]; p < rowPtr[i + 1]; ++p)
            sum += values[p] * x[colIdx[p]];
        y[i] = sum;
    }
}

int main() {
    /* 5 x 5 稀疏矩阵，7 个非零元，按行优先给出三元组 */
    const Triple coo[] = {
        {0, 0, 4}, {0, 2, 6},
        {1, 1, 5},
        {2, 0, 1}, {2, 3, 7},
        {3, 2, 2},
        {4, 4, 3}
    };
    const int nnz = 7, rows = 5;
    const int x[5] = {1, 2, 3, 4, 5};
    int y[5];

    cooToCsr(coo, nnz, rows);

    printf("values  :");
    for (int p = 0; p < nnz; ++p) printf(" %2d", values[p]);
    printf("\ncolIdx  :");
    for (int p = 0; p < nnz; ++p) printf(" %2d", colIdx[p]);
    printf("\nrowPtr  :");
    for (int i = 0; i < rows; ++i) printf(" %2d", rowPtr[i]);
    printf(" %2d", nnz);          /* 哨兵 rowPtr[rows] */
    printf("\n");

    spmvCsr(rows, x, y);
    printf("x       :");
    for (int i = 0; i < rows; ++i) printf(" %2d", x[i]);
    printf("\ny = A*x :");
    for (int i = 0; i < rows; ++i) printf(" %2d", y[i]);
    printf("\n");
    return 0;
}

/* 运行结果（g++ -std=c++17 -O2，与图 6-12 完全对应）：
     values  :  4  6  5  1  7  2  3
     colIdx  :  0  2  1  0  3  2  4
     rowPtr  :  0  2  3  5  6  7
     x       :  1  2  3  4  5
     y = A*x : 22 10 29  6 15
   手工核对第 0 行：y[0] = values[0]*x[colIdx[0]] + values[1]*x[colIdx[1]]
                       = 4*x[0] + 6*x[2] = 4 + 18 = 22  ✔
   注意 rowPtr[] 的长度是 rows + 1：多出来的那个哨兵 nnz = 7 让每行的区间
   都能写成半开区间，内层循环因此少一个 if —— 这就是 CSR 比三元组表简洁的地方。 */