简介面向 C 课程实验与面向对象编程练习者的矩阵类项目资料围绕 4×4、double 类型矩阵封装展开可用于理解类设计、构造函数、运算符重载与矩阵运算实现。内容包括默认构造为单位阵、拷贝构造、二维数组初始化以及加、减、乘、幂次、输入输出、赋值和下标访问等重载练习并实现求逆与转置适合课堂实验提交、课后复现和查漏补缺。压缩包共 7 个文件约 158KB涵盖 cpp 源码、sln 与 vcxproj 项目配置、docx 实验报告、vsdx 图示及 filters 等工程辅助文件便于直接用 Visual Studio 打开梳理结构。已有 139 人学习可作为矩阵运算与 C 运算符重载的完整参考案例。1. 一个 4×4 矩阵类能覆盖 C 的多个关键能力有人拿着自己写的 4×4 求逆来找我代码两百行编译零警告跑出来的 A·A⁻¹ 对角线是 1.0000003非对角线是 0.42。他以为是代数余子式抄错了实际上是消元时主元已经掉到 1e-17代码还在继续往下除。4×4 矩阵在 C 矩阵练手项目里算得上是性价比最高的一个问题初始化、转置、求逆、访问四个入口能把构造语义、const 正确性、运算符重载、浮点稳定性一次性串完。这份 C 实验项目适合刚学完类与运算符、想找一个真能跑出数的小工程的人也适合手头图形或物理模块里那个东拼西凑的 Matrix4 需要重捋一遍的老手。下面按我拆它的顺序来先骨架再核心算法然后怎么验最后怎么改成自己的东西。2. 工程骨架与类设计从 Matrix.h 到 CMake 的最小可编译闭环拿到一份矩阵类的代码我第一眼看的不是求逆那段而是头文件和目录。头文件里如果能看出「谁拥有数据、谁负责检查、谁负责抛异常」后面八成的 bug 就已经被挡住了。4×4 这个规模很微妙——小到可以完全展开写死大到如果不设计接口就会到处复制粘贴所以骨架怎么搭直接决定了这份代码是「交完作业就扔」还是「能放进项目里用」。2.1 目录结构与构建四个文件起步这个体量的项目不需要复杂分层但一定要把声明和实现分开否则测试代码没法独立编译。我一般按下面这个结构铺编译器和 IDE 都不会挑毛病matrix4x4/ ├── CMakeLists.txt ├── include/ │ └── Matrix.h # 类声明 内联小函数 ├── src/ │ └── Matrix.cpp # 构造、求逆、行列式等重逻辑 └── tests/ └── test_matrix.cpp # 不依赖第三方框架的断言测试对应的 CMakeLists.txt 只做三件事定标准、编库、连测试。cmake_minimum_required(VERSION 3.16) project(matrix4x4 LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) add_library(matrix STATIC src/Matrix.cpp) target_include_directories(matrix PUBLIC include) add_executable(matrix_test tests/test_matrix.cpp) target_link_libraries(matrix_test PRIVATE matrix) enable_testing() add_test(NAME matrix_test COMMAND matrix_test)这里几个参数值得说清楚。CMAKE_CXX_EXTENSIONS OFF关掉编译器私有扩展保证-stdc17是严格模式避免在某个编译器上能过、换个环境就报错。target_include_directories用PUBLIC而不是PRIVATE是因为测试目标要#include Matrix.h库的 include 路径得传递出去。add_library(... STATIC)打成静态库链接时不用管动态库路径本地跑测试最省事。编译跑测试三步cmake -S . -B build -DCMAKE_BUILD_TYPEDebug cmake --build build -j ctest --test-dir build --output-on-failureDebug 构建很关键。Release 下编译器可能把浮点中间结果放进 80 位寄存器某些误差用例的表现会和 Debug 不一致排查数值问题时先确认自己在哪个构建类型下。2.2 成员布局为什么用 std::arraydouble,16 行优先存储方式有三种常见写法double m[4][4]、std::arraystd::arraydouble,4,4、std::arraydouble,16。我选第三种理由不玄学第一内存连续且类型完整。double m[4][4]作为函数参数会退化成double(*)[4]sizeof直接失效std::arraydouble,16传参时不会退化能按值传递也能按引用传递还能直接memcpy给图形接口。第二行优先row-major跟大多数数学教材的下标习惯一致m[r*4c]一眼能对上A(r,c)如果后面要交给某些列优先的图形 API只在传数据那一步做一次转置不在类内部埋坑。第三16 个元素定长std::array没有堆分配构造和拷贝都是纯栈操作塞进结构体数组里也不会产生指针间接寻址。// include/Matrix.h #pragma once #include array #include cstddef #include initializer_list #include stdexcept class Matrix { public: static constexpr int N 4; Matrix() default; // 全零 explicit Matrix(double diag); // 对角阵 Matrix(std::initializer_listdouble init); // 行优先 16 个元素 static Matrix identity(); static Matrix zero(); double operator()(int r, int c); const double operator()(int r, int c) const; double* data() { return m_.data(); } const double* data() const { return m_.data(); } private: std::arraydouble, N * N m_{}; // 行优先默认全零 };std::arraydouble,16 m_{}后面这对花括号不能省。std::array是聚合类型成员声明处不加{}初始化默认构造出来的 16 个 double 是未定义值你会在operator*里读出随机数而且只在某些编译选项下复现——这类问题最难查。operator()而不是operator[]是因为矩阵要两个下标。C23 之前operator[]只能收一个参数硬要做成m[4][4]风格的访问就得返回行代理对象多一层间接寻址还容易悬垂。用operator()一次到位代价只是括号看起来不像数组习惯就好。2.3 三个构造入口与静态工厂把「初始化」这件事收口初始化的写法最容易发散。有人写Matrix m; m.set(0,0,1.0); ...十六行有人加一个void init(double a00, double a01, ...)十六个参数。这两种我都不推荐——前者啰嗦且容易漏后者参数顺序写反了编译器一声不吭。收口成三个入口就够了Matrix::Matrix(double diag) { for (int i 0; i N; i) m_[i * N i] diag; } Matrix::Matrix(std::initializer_listdouble init) { if (init.size() ! N * N) throw std::invalid_argument(Matrix: initializer_list must have 16 elements); std::copy(init.begin(), init.end(), m_.begin()); } Matrix Matrix::identity() { Matrix r; for (int i 0; i N; i) r.m_[i * N i] 1.0; return r; }Matrix(double diag)加了explicit防止Matrix m 2.0;这种隐式转换悄悄发生——标量乘矩阵是operator*(double)的活不该由构造函数代劳。initializer_list版本做长度校验少给或多给一个元素当场抛异常比后面算出一堆垃圾值强。identity()和zero()做成静态工厂调用点写Matrix::identity()比Matrix(1.0)语义清楚得多读代码的人不用回头翻构造函数。提示initializer_list构造函数不要和Matrix(double)同时存在歧义写法。Matrix{1.0}会优先匹配列表版本然后因为长度不对抛异常想构造对角阵必须写Matrix(1.0)圆括号。2.4 访问接口带检查的 operator() 和不带检查的 data()访问接口要同时满足两类调用方调试期的业务代码要安全性能敏感的循环要快。所以两个版本都给double Matrix::operator()(int r, int c) { if (r 0 || r N || c 0 || c N) throw std::out_of_range(Matrix: index out of range); return m_[r * N c]; } const double Matrix::operator()(int r, int c) const { if (r 0 || r N || c 0 || c N) throw std::out_of_range(Matrix: index out of range); return m_[r * N c]; }两个版本都要写这是 const 正确性的基本要求。只写非 const 版本const Matrix类型的引用和值就完全没法读元素operator*的const成员函数里调o(i,j)会直接编译失败。实现上让非 const 版本调用 const 版本再const_cast掉也可以四行代码不值得引入const_cast直接重复一遍更直白。data()是给底层用的逃生通道比如把结果直接喂给渲染 API 或者做memcpy序列化。给它就意味着放弃边界保护所以我在注释里写清楚「调用方自己保证 16 个元素的读写范围」。这种口子留一个就够不要再额外暴露at()、get()、value()三套等价接口API 面积越大越难维护。3. 核心运算实现求逆为什么必须走高斯-约当加部分选主元类骨架搭完真正吃时间的部分来了。转置和乘法是机械劳动行列式和求逆才是这份项目里唯一有算法含量的地方。很多人第一反应是伴随矩阵4×4 的逆等于伴随矩阵除以行列式公式在教材上白纸黑字写着照着展开 16 个代数余子式就行。这条路线能跑通但我不建议在工程里这么写原因和替代方案下面拆开讲。3.1 转置与乘法先定语义再写循环转置的语义要先定清楚是原地修改还是返回新对象。原地版签名是Matrix transposeInPlace()返回引用支持链式调用返回新版是Matrix transpose() constconst对象也能调。我的做法是两个都提供但把原地版命名成带InPlace后缀调用点一眼能看出这个对象被改了。Matrix Matrix::transpose() const { Matrix r; for (int i 0; i N; i) for (int j 0; j N; j) r.m_[j * N i] m_[i * N j]; return r; }循环顺序是「读原矩阵的 (i,j)写结果的 (j,i)」两次下标计算换位置就完成转置。4×4 规模下缓存不敏感怎么写都对但如果哪天扩成 NxN 模板写成按行扫描原矩阵、按列写结果访存局部性会好很多。乘法用「i-k-j」三层循环而不是教材上的「i-j-k」Matrix Matrix::operator*(const Matrix o) const { Matrix r; // 默认构造全零 for (int i 0; i N; i) { for (int k 0; k N; k) { const double aik m_[i * N k]; if (aik 0.0) continue; // 稀疏场景下跳过整行累加 for (int j 0; j N; j) r.m_[i * N j] aik * o.m_[k * N j]; } } return r; }i-k-j 顺序把内层循环变成对o的第 k 行和r的第 i 行的连续扫描aik提到最内层外面只读一次。if (aik 0.0) continue这一句在矩阵稀疏时能省掉整轮内层循环在稠密矩阵上也就是四次比较的开销不亏。参数o按 const 引用传避免 16 个 double 的拷贝。3.2 行列式LU 分解顺带把主元也算出来行列式如果用代数余子式展开4×4 要算 4 个 3×3 行列式每个再展开成 3 个 2×2代码量大且容易错一个符号。用带选主元的高斯消元二十行搞定顺便还能拿到「这个矩阵是否奇异」的信息double Matrix::determinant() const { double a[N][N]; for (int r 0; r N; r) for (int c 0; c N; c) a[r][c] m_[r * N c]; double det 1.0; int sign 1; for (int col 0; col N; col) { int pivot col; for (int r col 1; r N; r) if (std::abs(a[r][col]) std::abs(a[pivot][col])) pivot r; if (std::abs(a[pivot][col]) 1e-12) return 0.0; // 数值上视为奇异 if (pivot ! col) { for (int c 0; c N; c) std::swap(a[col][c], a[pivot][c]); sign -sign; // 交换两行行列式变号 } det * a[col][col]; for (int r col 1; r N; r) { const double f a[r][col] / a[col][col]; for (int c col; c N; c) a[r][c] - f * a[col][c]; } } return sign * det; }这里拷了一份a[N][N]出来因为消元是破坏性操作不能动原矩阵的m_。sign记录行交换次数每交换一次取反——这是行列式性质里最容易漏的一条漏了之后对角线元素的乘积是对的但符号随机错。内层for (int c col; ...)从col开始而不是 0因为第 col 列在消元后注定是 0不用算。阈值 1e-12 是经验值4×4 双精度下小于这个量级的主元基本等价于病态矩阵返回 0 比返回一个 1e-300 的数更有用。3.3 求逆高斯-约当消元加部分选主元的完整实现求逆的正统做法是在原矩阵右边拼一个单位阵凑成 4×8 的增广矩阵然后对左边 4 列做消元把左边化成单位阵右边就是逆Matrix Matrix::inverse() const { constexpr double EPS 1e-12; double a[N][2 * N]; // 拼增广矩阵 [A | I] for (int r 0; r N; r) { for (int c 0; c N; c) a[r][c] m_[r * N c]; for (int c 0; c N; c) a[r][N c] (r c) ? 1.0 : 0.0; } for (int col 0; col N; col) { // 1) 部分选主元在第 col 列剩下的行里找绝对值最大的 int pivot col; for (int r col 1; r N; r) if (std::abs(a[r][col]) std::abs(a[pivot][col])) pivot r; if (std::abs(a[pivot][col]) EPS) throw std::domain_error(Matrix::inverse: singular matrix); // 2) 换行 if (pivot ! col) for (int c 0; c 2 * N; c) std::swap(a[col][c], a[pivot][c]); // 3) 主元行归一化让 a[col][col] 1 const double inv 1.0 / a[col][col]; for (int c 0; c 2 * N; c) a[col][c] * inv; // 4) 消去其他所有行的第 col 列不只是下面几行 for (int r 0; r N; r) { if (r col) continue; const double f a[r][col]; if (f 0.0) continue; for (int c 0; c 2 * N; c) a[r][c] - f * a[col][c]; } } Matrix out; for (int r 0; r N; r) for (int c 0; c N; c) out.m_[r * N c] a[r][N c]; return out; }四步的先后顺序不能换。选主元必须在归一化之前否则除的就是一个接近零的数归一化必须在消元之前否则要额外处理除a[col][col]的除法。第 4 步的循环从r 0开始而不是col 1这是高斯-约当和高斯消元的区别——高斯消元只消下方得到上三角后还要回代约当形式把上方也一起消直接得到对角阵省掉回代那一步代价是循环次数翻倍。4×4 这个规模多出来的一百多次乘加完全可以忽略换来的是代码短、不容易写错。参数说明EPS是奇异性判定的阈值双精度下 1e-12 对 4×4 够用如果换成float存储要放宽到 1e-5 左右。抛std::domain_error而不是返回一个空矩阵或全零矩阵是因为「求逆失败」和「逆是零矩阵」在语义上完全不同用返回值区分会让调用方每次都要额外判断。异常信息里带上函数名出错时堆栈一眼能看到源头。注意返回的Matrix是值类型16 个 double 的拷贝在 4×4 规模下编译器基本会用 RVO 消掉。如果哪天扩成 100×100这里要改成bool inverse(Matrix out) const的写法把分配交给调用方。3.4 运算符重载与异常约定加减乘和标量乘按值返回参数 const 引用这些都是常识但有两处值得提前定规矩。第一不重载operator。浮点矩阵的相等判断永远有歧义用(A-B).norm() eps还是逐元素绝对误差取决于业务场景。我提供一个isApprox(const Matrix o, double eps)成员函数名字就把「这是近似比较」写在脸上了避免有人写出if (A B)然后在另一个编译器上结果飘。第二加减法不检查维度。4×4 定长类里维度恒定检查是多余的但如果后面要模板化成 NxN这条约定必须改成static_assert或者运行时抛异常否则 A(3×3) B(4×4) 会读越界。Matrix Matrix::operator(const Matrix o) const { Matrix r; for (int i 0; i N * N; i) r.m_[i] m_[i] o.m_[i]; return r; } Matrix Matrix::operator*(double s) const { Matrix r; for (int i 0; i N * N; i) r.m_[i] m_[i] * s; return r; }这两个函数直接在一维数组上做不用双重循环因为加减和标量乘本来就是逐元素的跟行列结构无关。写成i N*N一维遍历编译器更容易向量化。4. 测试与验证怎么证明求逆是对的矩阵类最大的陷阱是「看着对」。输出十六个数全是 0.9999 或者 1.0001 这种量级肉眼看不出错在哪。所以验证必须自动化而且要往数值稳定性上使劲不能只测一组数据。4.1 一个不引第三方库的断言宏实验项目没必要为了测试引入框架二十行就够#include cmath #include cstdio #include cstdlib #include algorithm static int g_failed 0; #define CHECK(cond) do { \ if (!(cond)) { \ std::fprintf(stderr, FAIL %s:%d %s\n, \ __FILE__, __LINE__, #cond); \ g_failed; \ } \ } while (0) static bool nearlyEqual(double a, double b, double eps 1e-9) { const double scale std::max({1.0, std::abs(a), std::abs(b)}); return std::abs(a - b) eps * scale; } static bool matricesClose(const Matrix x, const Matrix y, double eps 1e-9) { for (int i 0; i Matrix::N; i) for (int j 0; j Matrix::N; j) if (!nearlyEqual(x(i, j), y(i, j), eps)) return false; return true; }nearlyEqual用相对误差而不是绝对误差scale取max(1.0, |a|, |b|)保证在数值接近零时退化成绝对误差数值很大时是相对误差。这套判据跟常见的浮点比较写法等价不会出现「1e-10 和 2e-10 被判为不相等」这种过严的情况。4.2 数值验收A·A⁻¹ 与单位阵的误差怎么卡最核心的一条测试是「逆乘原矩阵应该等于单位阵」。误差阈值定多少有讲究4×4 双精度、元素量级在个位数时1e-12 是合理上限实际大多在 1e-15 附近。如果测出来是 1e-5那多半是算法路径有问题不是阈值太严。int main() { // 用例 1单位阵的逆是它自己 const Matrix I Matrix::identity(); CHECK(matricesClose(I.inverse(), I, 1e-15)); // 用例 2上三角阵手算可验 // A [[4,7,2,3],[0,5,1,2],[0,0,6,1],[0,0,0,8]]det 4*5*6*8 960 Matrix A{4,7,2,3, 0,5,1,2, 0,0,6,1, 0,0,0,8}; CHECK(nearlyEqual(A.determinant(), 960.0, 1e-12)); CHECK(matricesClose(A * A.inverse(), I, 1e-12)); CHECK(matricesClose(A.inverse() * A, I, 1e-12)); // 用例 3转置两次回到原点 CHECK(matricesClose(A.transpose().transpose(), A, 0.0)); // 用例 4det(A^-1) 应该等于 1/det(A) CHECK(nearlyEqual(A.inverse().determinant(), 1.0 / 960.0, 1e-12)); // 用例 5奇异矩阵必须抛异常 bool thrown false; try { Matrix S{1,2,3,4, 2,4,6,8, 0,0,0,0, 1,1,1,1}; // 第 2 行是第 1 行的 2 倍 (void)S.inverse(); } catch (const std::domain_error) { thrown true; } CHECK(thrown); if (g_failed 0) { std::puts(all tests passed); return 0; } std::fprintf(stderr, %d check(s) failed\n, g_failed); return 1; }用例 2 特意选了上三角阵行列式手算就能对上 960如果determinant()算出来是 -960 就说明sign逻辑写反了算出来是 4 或者 5 就说明内层循环的det *放错位置了。用例 5 的奇异矩阵故意让第二行是第一行的两倍这种整数构造的奇异阵在浮点下依然精确奇异测出来的结果稳定不像随机生成的近奇异矩阵那样受编译器影响。4.3 边界与病态用例清单跑通主流程只是开始下面这几类输入才是真正暴露问题的用例构造方式关注点单位阵Matrix::identity()逆是否原样返回det 是否为 1全零阵Matrix()必须抛异常不能除零对角阵Matrix(3.0)逆是对角线上的倒数快速验证近奇异阵主元取 1e-13是抛异常还是给出大数与 EPS 阈值是否一致大数量级阵对角线放 1e8归一化后是否溢出误差是否随条件数放大扰动矩阵单位阵加 1e-10 扰动逆的元素量级应该到 1e10验证没有误判奇异跑完这张表基本可以把「实现有 bug」和「数值本来就病态」区分开。前者表现为大片元素错得离谱后者表现为除了少数元素之外都准而那少数元素对应的是条件数很大的方向。判据是把误差跟条件数乘一下看是否落在机器精度范围内。5. 避坑与排查4×4 矩阵类最容易翻车的五个点下面五条有四个我亲眼见过别人踩还有一个是自己踩的。每条按「现象、原因、解决」写照着对号入座就行。5.1 现象const Matrix传进函数后编译报错现象一个只读的打印函数签名是void dump(const Matrix m)里面写m(0,0)直接报「no matching function for call to operator()」。原因只实现了非 const 版本的operator()const 对象只能调 const 成员函数编译器找不到匹配的重载。解决两个版本都实现const 版本返回const double。同样的坑还会出现在begin()、end()、data()这些访问器上凡是返回内部引用的接口都要配一份 const 版本。5.2 现象求逆结果每次运行都不一样Debug 和 Release 差得更远现象同一组输入两次运行的结果差在 1e-15 量级Release 下误差放大到 1e-6。原因一是std::arraydouble,16 m_;没写{}默认构造出来是未初始化内存构造函数只填了对角线其余位置是随机值二是 Release 下编译器做浮点重排-ffast-math或类似选项中间结果精度不同。解决成员声明处写std::arraydouble, N*N m_{};。排查是不是这个问题用Matrix m; m.data()[5]打出来看一眼如果是个奇怪的数那就是未初始化。数值相关的测试固定用 Debug 构建跑别在 Release 下对比精度。5.3 现象换个矩阵就除零崩溃或者结果出现 inf现象大部分矩阵求逆正常碰到某个特定矩阵inverse()里1.0 / a[col][col]直接得到 inf。原因没做行交换。上三角结构很好的矩阵不会出问题一旦第一行第一列恰好是 0比如交换矩阵[[0,1],[1,0]]扩展成 4×4主元为 0除法就炸了。解决老老实实做部分选主元在第 col 列剩余行里找绝对值最大的换上。这一步还有一个隐藏作用即使主元不为 0 只是很小选主元也能把它换成一个更合适的值显著降低误差放大。这是整份代码里唯一不能省的一步。5.4 现象A.inverse().inverse()断在中间或者结果莫名改变现象把逆矩阵写回原对象后续计算全乱。原因transpose()和inverse()都返回Matrix值类型临时对象的生命周期到分号结束。如果有人为了省事把它改成返回Matrix返回*this再配合auto r A.inverse();这种写法就变成了对临时对象取引用编译能过运行时读到已析构的内存。解决求逆、转置这类产生新值的操作一律按值返回不要返回引用。真要原地修改名字里带InPlace并且返回Matrix让调用点明确知道*this被改了。auto r A.inverse();是对的auto r A.inverse();应该被-Wreturn-stack-address之类的编译告警拦住把-Wall -Wextra打开。5.5 现象if (A B)编译不过或者偶尔为真偶尔为假现象想让两个矩阵比较相等发现要么没这个运算符要么加上之后结果不稳定。原因浮点矩阵的精确相等几乎永远为假A * A.inverse()的结果跟单位阵总有 1e-16 的差。重载operator会诱导调用方写出语义错误的比较。解决不重载提供bool isApprox(const Matrix o, double eps 1e-9) const用相对误差比较。函数名自带「近似」两个字调用方一看就知道自己在做什么。如果某些场景确实需要精确相等比如序列化后比对字节那就比memcmp(data(), o.data(), 16*sizeof(double))别混在数学比较里。6. 二次开发把 4×4 固长类改成 N×N 模板并接进渲染管线4×4 这个尺寸最适合拿来做模板化的第一步。把static constexpr int N 4提到模板参数上去改动量小收益却很直接3×3 的法线矩阵、2×2 的二维旋转、4×4 的变换都能用同一份实现而编译器会为每个 N 生成一份专门的代码循环边界全都在编译期确定。template int N class MatrixN { static_assert(N 1, MatrixN: N must be positive); public: MatrixN() default; explicit MatrixN(double diag) { for (int i 0; i N; i) m_[i * N i] diag; } MatrixN inverse() const; // 实现与 4x4 版逐行一致循环上界换成 N double determinant() const; MatrixN transpose() const; private: std::arraydouble, N * N m_{}; }; using Matrix4 MatrixN4; // 保留旧名字调用点零改动 using Matrix3 MatrixN3;迁移的诀窍是先加别名using Matrix4 MatrixN4;把现有代码全部切到别名上跑一遍测试确认行为一致再去删旧的Matrix类。上来就删类改调用点编译错误会多到没法判断哪一处是真的语义变了。接下来是性能。4×4 乘法在内层循环里固定跑 4 次编译器一般能展开但如果你要一帧算几千个骨骼矩阵值得写一个特化版本用double的四路累加代替内层循环template MatrixN4 MatrixN4::operator*(const MatrixN4 o) const { MatrixN4 r; for (int i 0; i 4; i) { for (int j 0; j 4; j) { // 四个乘加展开中间结果留在寄存器里 r.m_[i*4j] m_[i*40] * o.m_[0*4j] m_[i*41] * o.m_[1*4j] m_[i*42] * o.m_[2*4j] m_[i*43] * o.m_[3*4j]; } } return r; }这十六行是编译器最容易向量化的形态——内层是针对o同一列的四次连续读取-O2下通常能生成 SIMD 指令。改完必须重新跑第 4 章那套测试尤其是A * A.inverse()那条特化版本最容易在下标上写错一个位置然后determinant()是对的、乘法是错的非常具有迷惑性。验证改造是否成功我用的是一个笨办法但很有效随机生成一千组元素在 [-10, 10] 的矩阵两边都跑各自算A * A.inverse()与单位阵的最大偏差把最大值记下来。改造前后这个数应该在同一量级1e-15 到 1e-14如果改造后突然跳到 1e-8说明消元的循环上界或者下标哪里少算了一维。这个随机对拍脚本我建议留在仓库里每次动数学代码之前跑一遍比事后 debug 便宜太多。最后提醒一句关于「接进渲染管线」的现实问题很多图形接口和文件格式是列优先的而这份实现按行优先存。转换只需要在传数据的那一步做一次转置千万不要在类内部搞两套下标宏换来换去——我最早就是图省事加了#define AT(m,r,c) m[(c)*4(r)]结果一半函数用宏一半直接写调试了两天才发现某次新增的求逆函数漏了宏行列式在某个特定输入上总是差一个符号。从那以后我每次动这类数学代码都强制先跑随机对拍、再跑手算用例、最后才去看业务侧的现象——顺序反了就是在给自己找麻烦。希望帮到你。本文还有配套的精品资源点击获取
