简介：这份C++实验项目围绕4×4矩阵类Matrix_4x4的完整实现展开，面向正在学习C++面向对象、运算符重载与线性代数基础编程的高校学生及自学者，可用于课程实验、课后练习或项目复盘。压缩包共7个文件、约158KB，包含cpp源码、sln解决方案、vcxproj项目配置以及docx实验报告、vsdx图示等，便于在Visual Studio中直接打开、编译和对照理解工程结构。内容覆盖单位阵默认构造、拷贝构造、二维数组初始化、矩阵加减乘与幂次运算、输入输出、赋值和下标访问，并实现求逆与转置等核心功能，能帮助读者梳理类设计、重载规则和矩阵运算的实现思路。已有139人学习，适合作为矩阵类实验的参考样例与查漏补缺材料。

1. 一个 4x4 矩阵类，为什么值得单独拆一遍

做图形渲染的作业时，最玄学的 bug 往往不是着色器写错，而是模型矩阵的逆矩阵算出来带了一堆nan——法线变换一乘，整个模型直接黑掉或者变成马赛克。排查到最后，问题常常出在求逆函数没判断行列式接近 0 的情况，除零除出来的。这个 c++ 矩阵实验项目，做的就是把这件最基础也最容易翻车的事做扎实：封装一个 4x4 的Matrix类，覆盖单位矩阵与自定义数组两种初始化方式、矩阵加减乘、转置、行列式、求逆，以及operator()、operator[]这类元素访问接口。它适合刚学完类、构造析构、运算符重载，想找一个能真正编译跑起来、而不是只在 PPT 上看过语法的人；也适合需要回头补一遍线性代数代码手感的老手，把数值稳定性那几个坑重新踩一遍。

2. 矩阵类的骨架：从数据成员到接口设计

2.1 内存布局：为什么选double m[4][4]而不是double*

一个 4x4 的矩阵类，第一件事就是决定数据怎么存。可选方案有三种：double**的二级指针、std::vector<std::vector<double>>、以及定长的double m[4][4]。这个项目选了第三种，是个很务实的决定。

二级指针要自己管两次new，拷贝构造和析构都得多写一层循环，一旦漏掉就是内存泄漏；vector<vector<double>>虽然不用手动释放，但 4 个vector各自持有一块堆内存，行与行之间不连续，做矩阵乘的时候缓存命中率会明显打折。定长二维数组是连续存储的 16 个double，一共 128 字节，正好填满两条缓存行，编译器还能在开-O2之后把内层循环向量化。

// matrix4x4.h #pragma once #include <initializer_list> #include <stdexcept> class Matrix4x4 { public: static constexpr int N = 4; // 维度常量，避免代码里到处写魔数 4 Matrix4x4(); // 默认构造：单位矩阵 explicit Matrix4x4(double diag); // 对角构造：主对角线全为 diag Matrix4x4(const double (&arr)[N][N]); // 由二维数组拷贝构造 Matrix4x4(std::initializer_list<double> init); // 由 16 个花括号值构造 Matrix4x4(const Matrix4x4& other); // 拷贝构造 Matrix4x4& operator=(const Matrix4x4& other); // 赋值运算符 ~Matrix4x4() = default; // 无堆内存，默认析构即可 private: double m[N][N]{}; };

N定义成constexpr而不是宏，是为了后面扩展成模板时能直接拿来做非类型参数。析构函数写成= default是有意的信号：这个类没持有任何堆资源，拷贝就是纯粹的值语义，读代码的人一眼能看出不存在浅拷贝隐患。如果哪天有人把成员改成double* m，那这行= default就必须删掉并改成delete[]，这也是一个天然的改动提示点。

2.2 构造函数群：五种初始化入口的取舍

这个类的构造函数是整个设计里最值得看的部分，因为它把「矩阵怎么被创建出来」的几种典型场景都覆盖了。

#include "matrix4x4.h" // 默认构造：产出单位矩阵，避免用户拿到未初始化的垃圾值 Matrix4x4::Matrix4x4() { for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = (i == j) ? 1.0 : 0.0; } // 对角构造：Matrix4x4(2.0) 得到 2 倍单位阵，常用于缩放变换 Matrix4x4::Matrix4x4(double diag) { for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = (i == j) ? diag : 0.0; } // 数组构造：直接搬调用方给的 16 个数，行主序 Matrix4x4::Matrix4x4(const double (&arr)[N][N]) { for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = arr[i][j]; } // initializer_list 构造：方便在测试代码里写成 Matrix4x4{1,0,0,0, ...} Matrix4x4::Matrix4x4(std::initializer_list<double> init) { if (init.size() != N * N) throw std::invalid_argument("Matrix4x4 需要恰好 16 个元素"); auto it = init.begin(); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = *it++; } // 拷贝构造：定长数组直接逐元素复制，不需要 delete 旧值 Matrix4x4::Matrix4x4(const Matrix4x4& other) { for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = other.m[i][j]; } // 赋值运算符：先判自赋值，再逐元素覆盖 Matrix4x4& Matrix4x4::operator=(const Matrix4x4& other) { if (this == &other) return *this; // 自己赋值给自己，直接返回 for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) m[i][j] = other.m[i][j]; return *this; }

这里有几个容易被忽略的细节。第一，explicit加在double构造上，是为了阻止Matrix4x4 a = 3.0;这种隐式转换——虽然对数值类型看着无所谓，但一旦类被扩展成模板，隐式转换会带来一堆编译期歧义。第二，initializer_list版本里做了长度校验并抛std::invalid_argument，因为花括号初始化太容易被写错个数，宁可让它在构造时就炸掉，也不要留下一个半初始化的矩阵在后续计算里慢慢污染结果。第三，赋值运算符里的自赋值判断看着多余，但当矩阵类以后加了swap或移动语义，这个判断就是保命的。

成员初始化列表在这个类里用不上，因为数组没法在初始化列表里逐元素赋值。真正需要初始化列表的是把成员换成std::array<double, 16>的版本，那时候写: data{}就能全部清零，比在函数体里手写双重循环干净得多，这也是后面做模板化改造时的一个方向。

2.3 运算符重载骨架：加减乘与相等判断

矩阵类的手感，八成来自运算符重载写得好不好。这个项目实现了+、-、*、==、!=和复合赋值+=、-=、*=，全部是值语义返回新对象。

// 矩阵乘法：三重循环，i-k-j 顺序比 i-j-k 更缓存友好 Matrix4x4 operator*(const Matrix4x4& a, const Matrix4x4& b) { Matrix4x4 r(0.0); // 先构造全零矩阵 for (int i = 0; i < Matrix4x4::N; ++i) { for (int k = 0; k < Matrix4x4::N; ++k) { const double aik = a.get(i, k); // 把 a[i][k] 提到内层，减少重复取值 if (aik == 0.0) continue; // 稀疏场景下的快速跳过 for (int j = 0; j < Matrix4x4::N; ++j) r.set(i, j, r.get(i, j) + aik * b.get(k, j)); } } return r; } // 相等判断：浮点不能直接用 ==，用绝对误差阈值 bool operator==(const Matrix4x4& a, const Matrix4x4& b) { constexpr double eps = 1e-9; for (int i = 0; i < Matrix4x4::N; ++i) for (int j = 0; j < Matrix4x4::N; ++j) if (std::fabs(a.get(i, j) - b.get(i, j)) > eps) return false; return true; }

矩阵乘法把循环写成i-k-j而不是教科书上的i-j-k，是血泪经验：i-j-k在内层频繁跳行访问b，每次跨 32 字节，缓存几乎全失效；i-k-j让最内层连续读b的第k行，是顺序内存访问。数据规模只有 4x4 时差别肉眼看不出来，但这个习惯一旦养成，换成 64x64 或者 256x256 的矩阵就能差出好几倍。aik == 0.0的提前跳过是给稀疏矩阵留的口子，单位矩阵乘任何矩阵时能省掉四分之三的无效乘加，精度上没有副作用，因为乘 0 本来就是精确的。

==用误差阈值而不是裸==，是为了后面测A * A.inverse() == 单位矩阵时能通过。浮点求逆一定带舍入误差，如果坚持精确比较，测试会永远红着，最后逼得人去改测试而不是改实现，那就是本末倒置了。

3. 核心运算实现：求逆、转置与行列式的完整链路

3.1 元素访问：operator()、operator[]与边界检查

访问接口是这个类对外暴露最多的东西，写法直接决定了调用方代码好不好读。

// 非 const 版本：返回引用，可以出现在赋值号左边 double& Matrix4x4::operator()(int row, int col) { if (row < 0 || row >= N || col < 0 || col >= N) throw std::out_of_range("Matrix4x4 下标越界"); return m[row][col]; } // const 版本：返回 const 引用，供 const 对象和只读场景使用 const double& Matrix4x4::operator()(int row, int col) const { if (row < 0 || row >= N || col < 0 || col >= N) throw std::out_of_range("Matrix4x4 下标越界"); return m[row][col]; } // operator[] 只接受一个下标，按行主序映射到 m[row][col] double& Matrix4x4::operator[](int index) { if (index < 0 || index >= N * N) throw std::out_of_range("Matrix4x4 线性下标越界"); return m[index / N][index % N]; }

之所以同时提供两套，是因为使用场景不同。mat(1, 2)适合数学公式直译的代码，写mat(i, j)和线性代数教材上的 $a_{ij}$ 一一对应，可读性最好；mat[5]适合把矩阵当成 16 个连续浮点数来用的场景，比如往 GPU 传 uniform 变量、和第三方库的float[16]接口对接。两套接口最终都映射到同一块m内存，没有额外开销。

边界检查用抛异常而不是assert，这一点值得展开。assert在NDEBUG被定义的 Release 构建里会被整个裁掉，越界写入就变成静默的内存踩踏，症状会在几十行之外的地方冒出来，变成一个彻头彻尾的黑匣子。用throw std::out_of_range虽然带来一点点分支开销，但换来的是无论 Debug 还是 Release 都能定位到出错的那一行。真到了性能敏感的路径，再单独提供一套get(i, j)、set(i, j, v)的无检查接口，让调用方自己决定用哪个，这才是可控的取舍。

提示：operator[]的线性下标映射建立在行主序假设上。如果后面要和 OpenGL 的列主序矩阵对接，这层映射要么改成m[index % N][index / N]，要么在传参前整体转置一次，别在两个地方各改一半。

3.2 转置与行列式：求逆之前必须先算的两件事

转置是这一堆运算里最简单的一个，但它有两种写法，服务于不同场景。

// 写法一：返回新对象，原矩阵不变，适合函数式调用 Matrix4x4 Matrix4x4::transposed() const { Matrix4x4 r; for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) r.m[i][j] = m[j][i]; // 行列互换 return r; } // 写法二：原地转置，只交换上三角和对角线之外的元素 void Matrix4x4::transposeInPlace() { for (int i = 0; i < N; ++i) for (int j = i + 1; j < N; ++j) std::swap(m[i][j], m[j][i]); // 每对元素只交换一次 }

原地版本的内层循环从j = i + 1开始，是因为对角线上m[i][i]交换自己没意义，下三角的元素在前面的迭代里已经被换过了，再换一次等于白做。4x4 的规模下这点差异无所谓，但把「只遍历上三角」这个写法记住，处理更大矩阵时能省一半的交换操作。

行列式用代数余子式展开来算，4x4 展开成 4 个 3x3，每个 3x3 再展开成 3 个 2x2。

// 3x3 行列式，作为 4x4 展开的辅助函数 static double det3(const double s[3][3]) { return s[0][0] * (s[1][1] * s[2][2] - s[1][2] * s[2][1]) - s[0][1] * (s[1][0] * s[2][2] - s[1][2] * s[2][0]) + s[0][2] * (s[1][0] * s[2][1] - s[1][1] * s[2][0]); } double Matrix4x4::determinant() const { double total = 0.0; for (int j = 0; j < N; ++j) { // 沿第一行展开 if (m[0][j] == 0.0) continue; // 零元素跳过，省一次 det3 double sub[3][3]; for (int r = 1; r < N; ++r) { // 行从 1 开始，跳过被展开的第一行 int cc = 0; for (int c = 0; c < N; ++c) { if (c == j) continue; // 列跳过当前展开项所在的列 sub[r - 1][cc++] = m[r][c]; } } double sign = (j % 2 == 0) ? 1.0 : -1.0; // 代数余子式符号 (-1)^(0+j) total += sign * m[0][j] * det3(sub); } return total; }

这段代码的核心是「沿第一行展开」这个选择。理论上沿任意一行或一列展开结果都一样，但沿第一行写起来下标最简单，符号是(-1)^j。sub数组的填充用双指针r、cc手动控制行列索引，跳过了第 0 行和第j列，得到的就是余子式矩阵。m[0][j] == 0.0的跳过判断在稀疏矩阵上很有效，比如单位矩阵的行列式只需要算一次det3，其余三次全被跳过。

代数余子式展开的复杂度是 $O(n!)$，4x4 只需要 4 次 3x3 展开，完全够用。但如果哪天要扩展到 8x8 或者更大，这条路会瞬间爆炸，必须换成高斯消元法算行列式，复杂度降到 $O(n^3)$。这个类固定 4x4，所以选简单的写法没问题，认清边界在哪才是关键。

3.3 求逆的两套实现：伴随矩阵法与高斯-约当消元

求逆是这个项目里代码量最大、也最容易出错的部分，所以实现了两套：伴随矩阵法直白但慢，高斯-约当消元快且数值更稳，两套互相验证。

先看伴随矩阵法。核心公式是 $A^{-1} = \frac{1}{\det(A)} \cdot \text{adj}(A)$，其中伴随矩阵的每个元素是对应的代数余子式再转置。

Matrix4x4 Matrix4x4::inverseByAdjugate() const { double det = determinant(); if (std::fabs(det) < 1e-12) throw std::domain_error("矩阵奇异，不可逆"); // 关键：先判奇异 Matrix4x4 adj; for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { // 求 m[i][j] 的代数余子式：划掉第 i 行第 j 列后的 3x3 行列式 double sub[3][3]; int rr = 0; for (int r = 0; r < N; ++r) { if (r == i) continue; int cc = 0; for (int c = 0; c < N; ++c) { if (c == j) continue; sub[rr][cc++] = m[r][c]; } ++rr; } double sign = ((i + j) % 2 == 0) ? 1.0 : -1.0; // 注意这里是 adj[j][i] 而不是 adj[i][j]，转置体现在下标顺序上 adj.m[j][i] = sign * det3(sub); } } // 伴随矩阵整体除以行列式 for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) adj.m[i][j] /= det; return adj; }

这段里最容易写错的就是adj.m[j][i]。伴随矩阵的定义是「代数余子式矩阵的转置」，如果在填充时写成adj.m[i][j]，得到的是余子式矩阵本身，求出来的逆矩阵正好是正确结果的转置——这个 bug 在测试里表现为A * A.inverse()不等于单位矩阵，但等于单位矩阵的转置（其实单位矩阵转置还是自己，所以更阴险，对角矩阵测试根本发现不了，必须用非对称矩阵才能暴露）。

真正用在生产代码里的是高斯-约当消元。思路是把 $A$ 和单位矩阵拼成 $4 \times 8$ 的增广矩阵，对左边做行变换化成单位矩阵，右边自然就是逆矩阵。

Matrix4x4 Matrix4x4::inverse() const { double a[N][2 * N]; // 增广矩阵：左 A，右 I for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) a[i][j] = m[i][j]; for (int j = 0; j < N; ++j) a[i][j + N] = (i == j) ? 1.0 : 0.0; } for (int col = 0; col < N; ++col) { // 第一步：选主元，取当前列绝对值最大的行，减小舍入误差 int pivot = col; for (int r = col + 1; r < N; ++r) if (std::fabs(a[r][col]) > std::fabs(a[pivot][col])) pivot = r; if (std::fabs(a[pivot][col]) < 1e-12) throw std::domain_error("矩阵奇异，不可逆"); // 第二步：主元行换到当前行 if (pivot != col) for (int j = 0; j < 2 * N; ++j) std::swap(a[col][j], a[pivot][j]); // 第三步：当前行归一化，让主元变成 1 double d = a[col][col]; for (int j = 0; j < 2 * N; ++j) a[col][j] /= d; // 第四步：其他所有行消去当前列 for (int r = 0; r < N; ++r) { if (r == col) continue; double f = a[r][col]; if (f == 0.0) continue; for (int j = 0; j < 2 * N; ++j) a[r][j] -= f * a[col][j]; } } Matrix4x4 inv; for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) inv.m[i][j] = a[i][j + N]; return inv; }

pivot那一步是整段代码的灵魂。如果不选主元，直接按col行归一化，遇到主元位置是 0 的矩阵会直接除零；遇到主元很小的矩阵（比如 1e-15），归一化会把其他元素放大到 1e15 量级，后面的消元全被舍入误差淹没。选列主元之后，每一步除数的绝对值都是当前列最大的，误差被压到最小。

奇异判断的阈值1e-12也不是拍脑袋。对于double（约 15~16 位有效数字），如果主元绝对值比这个还小，说明矩阵条件数已经大到求出来的逆没有意义了。阈值设成1e-9会误伤一些本来可逆的病态矩阵，设成1e-15又挡不住真正的奇异矩阵，1e-12是这个规模下一个能用的经验值，实际项目里应该按矩阵的量级动态调整。

第四步用「消去所有行」而不是「只消去下方行」再回代，是因为消元法本身已经足够快，省掉回代那段代码换来的是逻辑更直白、更不容易写错。4x4 的规模下这点性能差别完全可以忽略，可维护性才是第一位的。

4. 编译、测试与运行：把 Matrix 类真正跑起来

4.1 目录结构与两种构建方式

这个实验项目的文件不多，但目录划分是标准的 C++ 工程结构，拿到手就能直接编译。

matrix4x4/ ├── include/ │ └── matrix4x4.h # 类声明，只放接口 ├── src/ │ ├── matrix4x4.cpp # 构造、访问、四则运算实现 │ ├── inverse.cpp # 行列式、转置、求逆实现 │ └── main.cpp # 演示入口 ├── tests/ │ └── test_matrix.cpp # 断言式测试 └── CMakeLists.txt

源文件按职责拆成两个.cpp，是有意的。matrix4x4.cpp里放的是没有数学风险的日常操作，inverse.cpp里放的是数值敏感的部分。将来要换一套求逆算法、或者接一个第三方线性代数库，只需要动一个文件，测试也能按文件粒度单独跑。

不用 CMake 的话，一行命令就能编出来，这是最低门槛的验证方式：

# 单文件直接编译，-O2 开优化，-Wall -Wextra 打开全部警告 g++ -std=c++17 -O2 -Wall -Wextra \ src/matrix4x4.cpp src/inverse.cpp src/main.cpp \ -Iinclude -o matrix_demo ./matrix_demo

-Wall -Wextra这两个参数千万别省。矩阵类里最容易出的一类低级错误是有符号和无符号比较（int和size_t混用），-Wextra会直接把它标出来；另一类是构造函数里成员没全部初始化，编译器也会警告。工程里习惯用 CMake 的话，配置长这样：

cmake_minimum_required(VERSION 3.16) project(matrix4x4 LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) # 调试构建开 -O0 -g，发布构建开 -O2，别在 Debug 里测性能 if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release) endif() add_library(matrix STATIC src/matrix4x4.cpp src/inverse.cpp) target_include_directories(matrix PUBLIC include) target_compile_options(matrix PRIVATE -Wall -Wextra) add_executable(matrix_demo src/main.cpp) target_link_libraries(matrix_demo PRIVATE matrix) enable_testing() add_executable(test_matrix tests/test_matrix.cpp) target_link_libraries(test_matrix PRIVATE matrix) add_test(NAME matrix_tests COMMAND test_matrix)

把核心实现编成STATIC库、可执行文件和测试都链接它，是为了避免每个可执行文件各自重新编译一遍源文件。target_include_directories用PUBLIC而不是PRIVATE，是因为测试和 demo 都要#include "matrix4x4.h"，头文件路径得往下传。

4.2 测试用例设计：从单位矩阵到随机矩阵

一个只有断言、不依赖第三方框架的测试文件，反而是最容易搬到别处的。

#include "matrix4x4.h" #include <cassert> #include <cmath> #include <cstdio> #include <random> static bool nearlyEqual(const Matrix4x4& a, const Matrix4x4& b, double eps = 1e-9) { for (int i = 0; i < Matrix4x4::N; ++i) for (int j = 0; j < Matrix4x4::N; ++j) if (std::fabs(a(i, j) - b(i, j)) > eps) return false; return true; } int main() { // 用例 1：单位矩阵求逆还是自己 Matrix4x4 I; assert(nearlyEqual(I.inverse(), I)); // 用例 2：对角矩阵求逆，对角线取倒数 Matrix4x4 D(2.0); Matrix4x4 expectedD(0.5); assert(nearlyEqual(D.inverse(), expectedD)); // 用例 3：非对称矩阵，专门用来抓「伴随矩阵忘了转置」的 bug Matrix4x4 A{ 4, 7, 2, 3, 1, 6, 5, 8, 9, 0, 3, 2, 5, 4, 1, 7 }; Matrix4x4 invA = A.inverse(); assert(nearlyEqual(A * invA, I)); // A 乘 A 逆应等于 I assert(nearlyEqual(invA * A, I)); // 交换顺序也要成立 // 用例 4：伴随矩阵法与高斯消元法结果互相印证 assert(nearlyEqual(A.inverseByAdjugate(), invA, 1e-7)); // 用例 5：奇异矩阵必须抛异常，而不是返回垃圾值 Matrix4x4 S{ 1, 2, 3, 4, 2, 4, 6, 8, // 第二行是第一行的两倍，行列式为 0 0, 1, 0, 0, 0, 0, 1, 0 }; bool thrown = false; try { S.inverse(); } catch (const std::domain_error&) { thrown = true; } assert(thrown); std::puts("all matrix tests passed"); return 0; }

用例 3 是整个测试里最关键的一条。如果只测单位矩阵和对角矩阵，前面说的「伴随矩阵忘记转置」的 bug 完全测不出来，因为这两类矩阵转置之后还是自己。必须用一个非对称、行列式不为 0 的矩阵，才能让A * invA和I的偏差暴露出来。

用例 4 做的是两套实现的交叉验证。两个算法独立写出来，在同一批输入上结果一致，才能说明至少没有出现系统性的逻辑错误。这里的误差阈值放宽到1e-7，因为伴随矩阵法经过了 16 次det3计算和一次大除法，累积误差本来就比高斯消元法大一个数量级。阈值不是越小越好，贴着算法本身的误差量级来设才合理。

用例 5 测的是异常路径。很多时候测试只覆盖「正常输入」，奇异矩阵、全零矩阵、含nan的矩阵这些边界情况完全没测，等上线以后在真实数据上撞见就直接崩。把异常路径写进测试，等于把这类问题的发现时间提前到了编译之后的第一秒。

4.3 数值校验：用随机矩阵做一轮批量验证

手写用例覆盖的是典型情况，真正要建立信心，还得靠随机矩阵批量跑。

int main() { std::mt19937 rng(42); // 固定种子，保证结果可复现 std::uniform_real_distribution<double> dist(-10.0, 10.0); int failed = 0; for (int t = 0; t < 1000; ++t) { double raw[4][4]; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) raw[i][j] = dist(rng); Matrix4x4 M(raw); if (std::fabs(M.determinant()) < 1e-6) continue; // 跳过近奇异的样本 Matrix4x4 inv = M.inverse(); // 计算 A * A^-1 - I 的最大元素误差 Matrix4x4 prod = M * inv; double maxErr = 0.0; for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) { double target = (i == j) ? 1.0 : 0.0; maxErr = std::max(maxErr, std::fabs(prod(i, j) - target)); } if (maxErr > 1e-8) { std::printf("case %d failed, maxErr = %e\n", t, maxErr); ++failed; } } std::printf("failed %d / 1000\n", failed); return 0; }

随机种子固定成 42，是为了让失败用例可以复现。随机测试最怕的就是「跑一次挂、再跑一次又好」，种子固定之后每次生成的矩阵完全一样，挂了就能直接定位到具体的t值，把那组数据打印出来单独调试。

跳过|det| < 1e-6的样本，是因为这类矩阵条件数很大，求出来的逆本来就带着可观的误差，拿它们去卡1e-8的精度阈值没有意义。真正的奇异判断在inverse()内部已经用1e-12挡过一道，这里只是把「接近奇异但没到阈值」的灰色地带排除在精度测试之外。

maxErr而不是布尔判断，是为了看误差的分布。如果 1000 个用例里有 900 个误差都在 1e-15 量级，只有几个在 1e-9，那说明算法没问题，是个别病态矩阵导致的；如果误差普遍在 1e-7 上下浮动，那就得怀疑主元选取或者归一化的顺序有问题。看数字比看红绿更能定位问题。

5. 避坑记录：4x4 矩阵类最常翻车的五个地方

5.1 求逆结果全是对的，但和正确矩阵差一个转置

现象：用对角矩阵、单位矩阵测试全部通过，换成一般的旋转矩阵或者随机矩阵，A * A.inverse()得到的不是单位矩阵，而是单位矩阵的转置——但单位矩阵转置还是自己，所以测试依然绿。真正在业务里用的时候，坐标变换的结果整体偏了一个角度。

原因：伴随矩阵法填充时写成了adj.m[i][j] = sign * det3(sub)，漏掉了伴随矩阵定义里那层转置。因为伴随矩阵是「代数余子式矩阵的转置」，正确写法应该是adj.m[j][i]。这个 bug 之所以阴险，是因为它只在对角矩阵上不可见，而大多数人随手写的第一个测试用例就是单位矩阵或对角矩阵。

解决：测试用例里必须包含至少一个非对称、非对角的矩阵，并且同时校验A * invA和invA * A。再进一步，用伴随矩阵法和高斯消元法互相印证，两套独立实现的结果对不上，就一定有一边写错了。

5.2 求逆函数在奇异矩阵上返回一堆nan或inf

现象：跑批量数据的时候，某个模型突然整个消失，或者渲染成一片黑色。打日志发现变换矩阵里全是nan，往上追溯是求逆函数返回了非法值，而且没有任何报错。

原因：高斯消元时直接拿a[col][col]做除数，没判断它是不是接近 0。奇异矩阵的主元位置恰好是 0，除下去就是inf；接近奇异的主元是一个极小值，除完之后其他元素被放大到 1e15 量级，再参与后续运算就溢出成nan。更麻烦的是nan有传染性，一次乘法就能污染整条计算链路。

解决：在归一化之前判fabs(a[pivot][col]) < 1e-12并抛std::domain_error。调用方要么在 try-catch 里处理，要么在调用前先用determinant()判一下。关键是别用「返回一个全零矩阵」这种方式糊弄过去，那只是把nan换成了更难排查的静默错误。

5.3 用assert做边界检查，Release 构建下越界写入没人管

现象：Debug 构建下跑得好好的，切到 Release 加-O2 -DNDEBUG之后，偶发崩溃，崩溃位置和矩阵代码隔着好几个函数，堆栈完全对不上。

原因：访问接口里写的是assert(row >= 0 && row < N)。NDEBUG一开，整个assert被预处理器抹掉，越界访问就变成了对栈上相邻内存的随机覆盖。4x4 矩阵只有 128 字节，越界写一格可能踩到旁边的局部变量，也可能踩到返回地址，症状完全随机。

解决：边界检查用throw std::out_of_range而不是assert。如果担心性能，就再提供一套get()、set()的无检查接口给热点路径，让调用方明确选择「我要快，我自己保证不越界」，而不是让编译器替所有人做决定。

5.4 矩阵乘法循环顺序写反，4x4 看不出问题，扩到 64x64 慢十倍

现象：小规模测试里性能正常，拿去做一批顶点变换，帧率明显掉下来，profile 一看热点全在矩阵乘法。

原因：循环写成for i / for j / for k，最内层每次访问b[k][j]都是跨行跳 32 字节，缓存行利用率接近零。4x4 的时候整个矩阵才 128 字节，全在 L1 里，怎么跳都快；矩阵一大，每次内层迭代都是一次缓存未命中。

解决：改成for i / for k / for j，最内层连续访问b的第k行，是顺序内存。顺手把a[i][k]提到内层循环外面存成局部变量，还能少一次寻址。这个习惯在小矩阵上看不出收益，但它是唯一一个不用改算法、只改循环顺序就能拿到数倍提升的地方。

5.5 拷贝构造和赋值运算符只写了一个，编译过了但行为不对

现象：Matrix4x4 b = a;正常，b = a;之后a的数据莫名其妙被改了，或者两个对象共享了同一块内存，改一个另一个跟着变。

原因：如果类里用了double*而不是定长数组，只写拷贝构造不写赋值运算符，编译器生成的默认赋值运算符做的是浅拷贝，两个对象指向同一块堆内存。析构的时候还会二次释放。这个项目用的是double m[4][4]，默认生成的拷贝构造和赋值本来就是逐元素复制，不会出问题；但一旦有人为了「支持动态维度」把成员改成指针，这个雷就埋下了。

解决：把「三件套」写全——拷贝构造、赋值运算符、析构。赋值运算符里先判自赋值if (this == &other) return *this;，再释放旧资源、分配新资源、复制数据。用定长数组的版本可以全部= default，但每次改成员类型时都要回头检查这三行还在不在。

6. 二次开发：把 4x4 写死改成模板维度

这个类最明显的局限就是维度写死在N = 4。想拿去做 3x3 的法线矩阵、2x2 的图像变换，就得复制一份代码改数字。改造方向很直接：把维度提成模板参数。

template <int N> class Matrix { public: Matrix() { // 默认构造为单位矩阵 for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) data_[i][j] = (i == j) ? 1.0 : 0.0; } double& operator()(int r, int c) { if (r < 0 || r >= N || c < 0 || c >= N) throw std::out_of_range("Matrix 下标越界"); return data_[r][c]; } template <int M> Matrix<M> resize() const { // 维度转换，多出的部分补单位阵 static_assert(M >= N, "目标维度不能小于源维度"); Matrix<M> r; for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) r(i, j) = data_[i][j]; return r; } private: double data_[N][N]{}; }; using Matrix4x4 = Matrix<4>; // 保留原来的名字，老代码不用改 using Matrix3x3 = Matrix<3>;

using别名是这次改造里最省事的一步：原来的Matrix4x4名字保持不变，所有调用点一行都不用动，编译期自动绑到Matrix<4>上。static_assert挡住的是Matrix<4>().resize<2>()这种把大矩阵往小矩阵塞的调用，让它直接编译失败而不是运行时截断。

改造过程中有几个验证点值得固定下来。第一，模板化之后必须把原来那 1000 组随机矩阵测试原封不动再跑一遍，确认Matrix<4>的行为和原来的Matrix4x4完全一致，这是改造没引入回归的唯一证据。第二，测一次Matrix<1>和Matrix<2>这种极小维度，因为模板代码里很多循环边界在N = 1时会退化，比如原地转置的内层循环j = i + 1直接不执行，det3那种辅助函数则完全用不上。第三，如果后面要加移动构造，记得把赋值运算符里的自赋值判断保留下来，模板化之后对象可能被std::move到容器里，自赋值的概率比想象中高。

求逆那一行的改动最大。伴随矩阵法里的det3是专门给 4x4 展开写的，模板化之后必须换成通用实现——要么递归展开（复杂度爆炸，只适合极小维度），要么干脆只用高斯-约当消元。我一般会保留伴随矩阵法但只在N == 4时通过if constexpr启用，其余维度全部走消元法，测试里两个都跑，对不上就报错。

最后说个教训。早些年写这类矩阵类，我总觉得「4x4 而已，随便写写就行」，直到有一次批处理工具里一个求逆返回了nan，污染了整批数据，而那个nan在流程里传了七八个模块才被发现。从那以后我每次动数值代码，都会先在测试里加一条A * A.inverse() == I的断言，再把奇异矩阵的异常路径跑一遍，才敢往下写业务逻辑。数值代码的坑不会在写的时候疼，只会在用的时候疼，提前花十分钟把这两条断言补上，能省掉后面一整天翻日志。希望帮到你。

本文还有配套的精品资源，点击获取