M01 浮点数与误差数值计算的世界观在写任何算法之前先理解这台机器的数值世界观。本模块地图数学原理部分建立 IEEE 754 与条件数两块基石回答「误差从哪来、问题本身有多敏感」算法设计部分给出评判算法的金标准——后向稳定性以及灾难性抵消的稳定化改写工具箱库级实现把这一切写进 matlib 的第一个模块并上机实测工业级解析看 Eigen 与 NumPy 如何暴露同样的机器参数以及为什么log1p、hypot早在 C99 就进了标准库。学习目标能说出 IEEE 754 double 的位结构并从位级解释0.1 0.2 ! 0.3能写出程序实测本机机器精度并与理论值 2⁻⁵² 对上能解释为什么「后向稳定 条件数」取代「前向误差小」成为评判算法的金标准能识别至少四种灾难性抵消场景并给出稳定化改写求和、二次求根、log1p、hypot能用 ULP 距离判断浮点相等说清比较浮点数为什么不可靠。1. 数学原理1.1 实数轴上有洞有限位宽装不下实数double 只有 64 位最多表示 2⁶⁴ ≈ 1.8×10¹⁹ 个不同的值而实数是无穷的。结论很简单绝大多数实数在机器里根本不存在每个浮点运算都是在「把真实结果搬到一个最近的可表示数」。对工程师意味着什么浮点误差不是 bug是这套表示系统的固有契约。我们要做的不是消灭误差而是知道误差有多大、什么时候会被放大。1.2 IEEE 754 double 的取值公式double 的 64 位分成三段段位数含义符号 s10 正 / 1 负偏移指数 e11实际指数 e − 10231023 为 bias尾数 f52小数部分规格化数隐含前导 1规格化数的取值公式(−1)s×2e−1023×1.f(-1)^s \times 2^{e-1023} \times 1.f(−1)s×2e−1023×1.f两个特殊区域要记住非规格化数e 0值为 (−1)^s × 2⁻¹⁰²² × 0.f填补 0 附近的空隙让下溢「逐渐发生」gradual underflow而不是断崖归零特殊值e 0x7FFf 0 时为 ±∞f ≠ 0 时为 NaN。例0.1 的二进制是无限循环小数 0.0001100110011…₂52 位尾数只能截断存储实际存下的是0.110⟶0.1000000000000000055511151231257827…0.1_{10} \longrightarrow 0.1000000000000000055511151231257827\ldots0.110⟶0.1000000000000000055511151231257827…这不是十进制转二进制的 bug而是「有限位宽装不下」的直接后果。1.3 ULP浮点数的间距随量级变化ULPUnit in the Last Place是相邻两个可表示浮点数之间的距离。在区间 [2^k, 2^{k1}) 内所有数的尾数部分均匀分布因此间距固定为ULP2k×2−522k−52\mathrm{ULP} 2^{k} \times 2^{-52} 2^{k-52}ULP2k×2−522k−52即指数每加 1间距翻倍1.0 附近的 ULP 是 2⁻⁵²10¹⁶ 附近的 ULP 已经是 2。对工程师意味着什么同一个绝对误差对小数值可能是致命的对大数值可能无关紧要。浮点世界必须用相对误差说话——这正是下一节的铺垫。1.4 机器精度 εmach 与舍入单元 u机器精度 εmach定义为 1.0 到下一个可表示数的距离。由上一节1.0 落在 [2⁰, 2¹) 区间所以εmach2−52≈2.22×10−16(double)\varepsilon_{\mathrm{mach}} 2^{-52} \approx 2.22 \times 10^{-16}\quad (\text{double})εmach2−52≈2.22×10−16(double)IEEE 754 默认的舍入到最近偶数模式下单次运算的相对误差不超过半个 ULP于是得到舍入单元uεmach22−53≈1.11×10−16u \frac{\varepsilon_{\mathrm{mach}}}{2} 2^{-53} \approx 1.11 \times 10^{-16}u2εmach2−53≈1.11×10−16推导任取 x ∈ [2^k, 2^{k1})半个 ULP 为 2^{k-53}相对误差最大出现在 x 取区间左端点时即 2{k-53}/2k 2⁻⁵³。对工程师意味着什么double 「免费」给你约 15–16 位十进制有效数字一次运算的损失不超过 u。后续所有算法的误差分析都以 u 为基本单位——这是本课程的「计量衡」。1.5 条件数问题自身的体质把「问题」看成函数 f输入 → 输出。输入的相对误差经过 f 会被放大多少倍这就是条件数κ(f;x)limδ→0∣δf/f(x)∣∣δx/x∣∣x f′(x)f(x)∣\kappa(f; x) \lim_{\delta \to 0} \frac{|\delta f / f(x)|}{|\delta x / x|} \left| \frac{x\, f(x)}{f(x)} \right|κ(f;x)δ→0lim∣δx/x∣∣δf/f(x)∣f(x)xf′(x)两个例子f(x) x²κ |x·2x/x²| 2。输入误差放大 2 倍良态f(x) 1 − xκ |x/(1−x)|。当 x 0.9999 时 κ ≈ 10⁴——输入的第 4 位有效数字会被放大成输出的第 1 位病态。对工程师意味着什么条件数是问题的先天体质与算法无关。病态问题上不存在「精确算法」——输入里根本没有那些信息。正确的姿势是先判断问题良态与否再谈算法算出大误差时先分清是「问题病态」还是「算法不佳」。2. 算法设计2.1 金标准后向稳定性衡量一个数值算法有两条路线前向误差算出的答案 f̂(x) 离真值 f(x) 有多远后向误差算出的答案是不是「某个 nearby 问题的精确解」——即是否存在 x̃ 使得 f̂(x) f(x̃)且 x̃ 离 x 很近。定义算法是后向稳定的当且仅当对每个输入 x都存在 x̃ x(1 O(u)) 使得算法输出恰为 f(x̃)。直观说法它「精确地解决了一个几乎一样的问题」。为什么后向稳定能成为金标准因为它和条件数组合后给出完整的前向误差界前向相对误差≲κ(x)⋅u\text{前向相对误差} \lesssim \kappa(x) \cdot u前向相对误差≲κ(x)⋅u即后向稳定算法 × 良态问题 可靠结果。这个框架还有一个工程上极好的性质——归因清晰结果不准要么问题病态κ 大换问题建模要么算法不后向稳定换算法。为什么不直接用「前向误差小」做标准前向误差同时由问题和算法决定库作者无法对任意用户输入做出承诺而后向稳定性是算法自身的性质一次证明、处处成立。这就是工业界选择它的原因——不是它更「强」而是它可认证、可归因。2.2 灾难性抵消减法不背锅数据先脏了两个相近的浮点数相减有效数字大量丢失这就是灾难性抵消catastrophic cancellation。定量地设 x̂ x(1δ₁)、ŷ y(1δ₂)|δ| ≤ u则∣误差∣∣x^−y^∣≤u⋅∣x∣∣y∣∣x−y∣\frac{|\text{误差}|}{|\hat x - \hat y|} \le u \cdot \frac{|x| |y|}{|x - y|}∣x^−y^∣∣误差∣≤u⋅∣x−y∣∣x∣∣y∣当 x ≈ y 时放大因子 ≈ 2|x|/|x−y|可以任意大。一个关键澄清减法本身是精确的。Sterbenz 引理保证当 x/2 ≤ y ≤ 2x 时x − y 没有舍入误差。抵消之所以「灾难」是因为相减的两个数在更早的运算里已经被污染减法只是把污染暴露出来。常见场景与稳定化改写本模块实现其中三个场景抵消来源稳定化改写二次方程求根 −b±√Δb² ≫ 4ac 时 √Δ ≈ |b|q 公式 韦达定理见 2.3方差 E[x²] − (E[x])²大均值小方差时两项相近Welford 在线算法课后练习(1x) − 1x 很小1 吞掉 x 的低位log1p(x)/expm1(x)√(x1) − √xx 大时两值相近有理化1/(√(x1)√x)(f(xh)−f(x))/hh 小时分子相近复步微分/符号微分毕业项目方向 3 的伏笔二次求根的稳定公式教科书公式 x (−b ± √Δ)/(2a) 在 b 与 √Δ 同号时必然抵消。稳定做法分两步构造 q −0.5·(b sign(b)·√Δ)保证括号内同号相加、绝不相消得一个远离抵消的根 x₁ q/a另一个根用韦达定理 x₁x₂ c/a 求出x₂ c/q同样避开减法。复杂度仍为 O(1)一次 sqrt、常数次四则运算代价几乎为零误差从「灾难级」降到 O(u)——M01 实验会实测这一差距。2.3 求和算法从 O(nu) 到 O(u)n 个数朴素累加的最坏误差界约为∣误差∣≲(n−1) u∑i∣xi∣|\text{误差}| \lesssim (n-1)\, u \sum_i |x_i|∣误差∣≲(n−1)ui∑∣xi∣误差随 n线性增长且与数据顺序有关。两个改进方向Kahan / Neumaier 补偿求和。核心思想用补偿变量 c 收集每一步被舍入丢掉的部分最后加回。经典 Kahan 假设「当前和总比新加入的数大」当数据顺序打破这一假设时会失效Neumaier 改进版补上了这个角落matlib 采用的就是它。伪代码s ← 0; c ← 0 for x in data: t ← s x if |s| ≥ |x|: c ← c ((s − t) x) # x 的低位被丢进这一项 else: c ← c ((x − t) s) # s 的低位被丢进这一项 s ← t return s c复杂度 O(n)每步多常数次浮点运算误差界降到约 2u·Σ|x|与 n 基本无关。备选方案对比为什么不选它们按绝对值排序后从小到大累加直觉正确但需要 O(n log n) 排序、必须拿到全部数据不适合流式且最坏误差界仍与 n 有关。弃用但思想保留pairwise两两分组求和把序列递归二分成树形求和误差 O(log n · u)无需维护补偿变量、缓存友好、可向量化。这是NumPy 的工业选择——它是「足够好 足够快」的折中细节见第 4 节与 M12。matlib 的选择流式精确场景用 Neumaierkahan_sum大规模向量化场景留待 M09 之后再考虑 pairwise。3. 库级实现3.1 matlib 骨架本模块为 matlib 建立目录与构建骨架code/ ├── CMakeLists.txt # C17MSVC 下开启 /utf-8 /W4 ├── matlib/ │ └── float_info.hpp # M01 数值底座header-only └── test_m01.cpp # 数值验证程序22 项断言约定从本模块开始生效C17、零第三方依赖、命名向 Eigen 看齐camelCase。3.2 float_info.hpp 关键接口位级视图。bits_of用 memcpy 做标准允许的位级双关decompose按 1|11|52 拆出三段inlinestd::uint64_tbits_of(doublex){std::uint64_tb;std::memcpy(b,x,sizeof(b));// memcpy 是标准允许的位级双关写法returnb;}机器精度实测。循环对半逼近直到 1 eps/2 被舍入回 1templatetypenameTTmachine_epsilon(){T epsT(1);for(intguard0;guard256;guard){// 防御式上限if(T(1)eps/2T(1))break;eps/2;}returneps;}ULP 比较。先把位模式映射成保序整数负数区间翻转、−0 归一到 0再做无符号差——用 unsigned 是为了覆盖 ±DBL_MAX 之间距离超过 int64 范围的极端边界inlinestd::int64_tordered_bits(doublex){conststd::int64_tbstatic_caststd::int64_t(bits_of(x));return(b0)?(INT64_MIN-b):b;}inlineboolnearly_equal(doublea,doubleb,std::uint64_tmax_ulp4);边界探测。从 1.0 出发倍增直到溢出、倍减直到即将进入非规格化区实测浮点世界的两堵墙2^1023 与 2⁻¹⁰²²实现见test_m01.cpp第 [8] 段。Neumaier 补偿求和kahan_sum、稳定求根solve_quadraticq 公式 韦达定理处理 a0 退化、Δ0、二重零根三个边界、稳定 hypothypot_stable按最大值缩放。完整代码见code/matlib/float_info.hpp。3.3 数值验证实测结果测试环境Windows x64MSVC 14.44VS 2022/std:c17 /O2 /EHsc /utf-8。完整验证程序为test_m01.cpp22 项断言全部通过。摘录关键输出eps_mach(double) 2.2204460492503131e-16 # 与理论值 2^-52 一致 eps_mach(float) 1.1920929e-07 # 2^-23 bits(0.1) 0x3FB999999999999A # 无限循环二进制被截断的证据 0.10.2 0.30000000000000004 0.3 0.29999999999999999, ulp 距离 1 # 只差 1 个 ULP naive_sum 0, kahan_sum 1000 # {1e16, 1×1000, -1e16}精确和 1000 sum(0.1 x 1e6): naive |err| 1.33e-06, kahan |err| 0 x²1e8x1 的小根 参考值 -1.0000000000000002e-08 naive 求根 -7.4505805969238281e-09 相对误差 0.255 # 错了 25% 稳定版求根 -1e-08 相对误差 1.65e-16 # 回到 u 量级 sqrt(3e200²4e200²) inf溢出 hypot_stable 5e200 ((1x)-1)/x 0 log1p(x)/x 1 (x1e-16) 边界探测: 最大 2 的幂 2^1023最小正规格化数 2^-1022 2^-1022 之下还有 52 级非规格化数直到 2^-1074 才归零——gradual underflow三个值得盯着看的结果朴素求和得 0ulp(10¹⁶) 2每个 1 恰好落在舍入的中点上被丢弃——不是「误差小」是「信息整体蒸发」kahan 求和误差为 010⁶ × fl(0.1) 与 100000 的差距约 5.55×10⁻¹²小于该量级的半个 ULP约 7.3×10⁻¹²最后的加法舍入恰好把它吸收回 100000.0naive 求根错 25%教科书公式在 b 10⁸ 时小根已完全不可信——这是「算法不稳定」与「问题病态」的区分样本问题求这个根本身良态是算法把误差放了进来。4. 工业级解析4.1 EigenNumTraits 暴露机器参数Eigen 把本模块手测的所有机器参数收进一个 traits 类NumTraitsT源码坐标Eigen/src/Core/NumTraits.hEigen 3.4.0epsilon()机器精度对应 matlib 的machine_epsilondummy_precision()经验容差epsilon 的小倍数用于近似比较——对应nearly_equal的思路但 Eigen 选了「容差倍数」路线我们选了「ULP 距离」路线后者在跨量级比较时更稳highest()/lowest()取值范围边界。后续模块中 Eigen 的各类阈值判断比如秩判定都从 NumTraits 取数——机器参数集中管理是库级设计的基本功。4.2 NumPyfinfo 与 pairwise 求和NumPy 侧的对应物NumPy 2.xnp.finfo(float)返回 eps / tiny / max分别对应机器精度、最小规格化数、最大可表示数np.spacing(1.0)返回 1.0 处的 ULP恰为 2⁻⁵²np.nextafter取相邻浮点数——正是 3.2 节ulp_distance的原料np.add.reduce即np.sum的底层不用朴素累加而是pairwise 求和NumPy 1.9 起源自社区对求和精度的长期讨论不超过 128 个元素的块内做展开直加超过则递归二分。实现位于源码numpy/_core/src/umath/loops.c.srcNumPy 2.0 之前为numpy/core/src/umath/loops.c.src可在源码中 greppairwise定位M12 会逐行拆。对照 Python 生态的另外两个选择math.fsum维护全部部分和精确但慢CPython 3.12 起内置sum()改用 Neumaier 补偿——与 matlib 的选择同款。4.3 C 标准库稳定化改写早就写进了标准C99 在math.h中加入log1p、expm1、hypot、fma——这些函数存在的唯一理由就是本模块第 2 节的稳定化改写。工业共识稳定形式是默认形式朴素写法是调用者的责任。4.4 差距分析教学版实现与工业版的差距不在算法思想本模块的公式与工业界完全同源而在工程维度无 SIMD、无多精度回退、无全格式覆盖float16/bfloat16、无平台特殊值穷举测试。这些维度分别对应后续模块SIMD → M09格式与 dispatch → M10–M12。5. 随堂实验目的在自己的机器上复现本模块全部结论——机器精度不是背下来的 2⁻⁵²是测出来的。验收数据要求实测 εmachfloat/double并与理论值对照构造一组数据使朴素求和误差可见而补偿求和精确复现二次求根的 naive vs 稳定版误差差给出hypot溢出边界的实测值。结论必须回答为什么 kahan 求和在「0.1 × 10⁶」实验里误差恰好为 0提示想清楚最终那次加法的舍入。详细步骤见lab.md将随「生成 M01 实验」指令产出。6. 本模块小结double 只有 2⁶⁴ 个值误差是表示系统的固有契约不是 bugULP 随指数翻倍浮点世界必须用相对误差说话εmach 2⁻⁵²舍入单元 u 2⁻⁵³是后续所有误差分析的计量单位条件数是问题的先天体质先判良态病态再谈算法金标准是后向稳定性前向误差 ≲ κ·u归因清晰、可认证灾难性抵消中减法本身精确Sterbenz锅在上游已被污染的数据Neumaier 补偿求和把求和误差从 O(nu) 压到 O(u)NumPy 用 pairwise 走了另一条同样成立的工业路线。在主线上的位置本模块四个阶段全部走通但重心在「数学原理 算法设计」——它提供的是贯穿全课程的世界观与误差语言库级实现刚打下 matlib 的地基真正的性能与架构主题从 M08 开始。7. 课后练习概念复盘εmach 2⁻⁵² 而相对舍入误差界 u 2⁻⁵³差一个因子 2。用自己的话解释这两个数各自的定义以及为什么「测 εmach」和「用 u 做误差分析」都合理。实现给 matlib 增加dot_stable用 Neumaier 补偿的向量点积。构造两组量级悬殊的向量实测 naive 点积与dot_stable的误差差。实现实现 Welford 在线方差算法单遍、数值稳定与「E[x²] − (E[x])²」两遍公式在「均值 10⁹、标准差 1」的数据上对比复现后者算出负方差的翻车现场。挑战选做找到 NumPyloops.c.src中的 pairwise 求和实现解释块阈值为什么取 128提示想想缓存行与展开因子并复现「pairwise 误差介于 naive 与 Kahan 之间」的曲线。 M01 ｜ 浮点数与误差：数值计算的世界观
M01 浮点数与误差数值计算的世界观在写任何算法之前先理解这台机器的数值世界观。本模块地图数学原理部分建立 IEEE 754 与条件数两块基石回答「误差从哪来、问题本身有多敏感」算法设计部分给出评判算法的金标准——后向稳定性以及灾难性抵消的稳定化改写工具箱库级实现把这一切写进 matlib 的第一个模块并上机实测工业级解析看 Eigen 与 NumPy 如何暴露同样的机器参数以及为什么log1p、hypot早在 C99 就进了标准库。学习目标能说出 IEEE 754 double 的位结构并从位级解释0.1 0.2 ! 0.3能写出程序实测本机机器精度并与理论值 2⁻⁵² 对上能解释为什么「后向稳定 条件数」取代「前向误差小」成为评判算法的金标准能识别至少四种灾难性抵消场景并给出稳定化改写求和、二次求根、log1p、hypot能用 ULP 距离判断浮点相等说清比较浮点数为什么不可靠。1. 数学原理1.1 实数轴上有洞有限位宽装不下实数double 只有 64 位最多表示 2⁶⁴ ≈ 1.8×10¹⁹ 个不同的值而实数是无穷的。结论很简单绝大多数实数在机器里根本不存在每个浮点运算都是在「把真实结果搬到一个最近的可表示数」。对工程师意味着什么浮点误差不是 bug是这套表示系统的固有契约。我们要做的不是消灭误差而是知道误差有多大、什么时候会被放大。1.2 IEEE 754 double 的取值公式double 的 64 位分成三段段位数含义符号 s10 正 / 1 负偏移指数 e11实际指数 e − 10231023 为 bias尾数 f52小数部分规格化数隐含前导 1规格化数的取值公式(−1)s×2e−1023×1.f(-1)^s \times 2^{e-1023} \times 1.f(−1)s×2e−1023×1.f两个特殊区域要记住非规格化数e 0值为 (−1)^s × 2⁻¹⁰²² × 0.f填补 0 附近的空隙让下溢「逐渐发生」gradual underflow而不是断崖归零特殊值e 0x7FFf 0 时为 ±∞f ≠ 0 时为 NaN。例0.1 的二进制是无限循环小数 0.0001100110011…₂52 位尾数只能截断存储实际存下的是0.110⟶0.1000000000000000055511151231257827…0.1_{10} \longrightarrow 0.1000000000000000055511151231257827\ldots0.110⟶0.1000000000000000055511151231257827…这不是十进制转二进制的 bug而是「有限位宽装不下」的直接后果。1.3 ULP浮点数的间距随量级变化ULPUnit in the Last Place是相邻两个可表示浮点数之间的距离。在区间 [2^k, 2^{k1}) 内所有数的尾数部分均匀分布因此间距固定为ULP2k×2−522k−52\mathrm{ULP} 2^{k} \times 2^{-52} 2^{k-52}ULP2k×2−522k−52即指数每加 1间距翻倍1.0 附近的 ULP 是 2⁻⁵²10¹⁶ 附近的 ULP 已经是 2。对工程师意味着什么同一个绝对误差对小数值可能是致命的对大数值可能无关紧要。浮点世界必须用相对误差说话——这正是下一节的铺垫。1.4 机器精度 εmach 与舍入单元 u机器精度 εmach定义为 1.0 到下一个可表示数的距离。由上一节1.0 落在 [2⁰, 2¹) 区间所以εmach2−52≈2.22×10−16(double)\varepsilon_{\mathrm{mach}} 2^{-52} \approx 2.22 \times 10^{-16}\quad (\text{double})εmach2−52≈2.22×10−16(double)IEEE 754 默认的舍入到最近偶数模式下单次运算的相对误差不超过半个 ULP于是得到舍入单元uεmach22−53≈1.11×10−16u \frac{\varepsilon_{\mathrm{mach}}}{2} 2^{-53} \approx 1.11 \times 10^{-16}u2εmach2−53≈1.11×10−16推导任取 x ∈ [2^k, 2^{k1})半个 ULP 为 2^{k-53}相对误差最大出现在 x 取区间左端点时即 2{k-53}/2k 2⁻⁵³。对工程师意味着什么double 「免费」给你约 15–16 位十进制有效数字一次运算的损失不超过 u。后续所有算法的误差分析都以 u 为基本单位——这是本课程的「计量衡」。1.5 条件数问题自身的体质把「问题」看成函数 f输入 → 输出。输入的相对误差经过 f 会被放大多少倍这就是条件数κ(f;x)limδ→0∣δf/f(x)∣∣δx/x∣∣x f′(x)f(x)∣\kappa(f; x) \lim_{\delta \to 0} \frac{|\delta f / f(x)|}{|\delta x / x|} \left| \frac{x\, f(x)}{f(x)} \right|κ(f;x)δ→0lim∣δx/x∣∣δf/f(x)∣f(x)xf′(x)两个例子f(x) x²κ |x·2x/x²| 2。输入误差放大 2 倍良态f(x) 1 − xκ |x/(1−x)|。当 x 0.9999 时 κ ≈ 10⁴——输入的第 4 位有效数字会被放大成输出的第 1 位病态。对工程师意味着什么条件数是问题的先天体质与算法无关。病态问题上不存在「精确算法」——输入里根本没有那些信息。正确的姿势是先判断问题良态与否再谈算法算出大误差时先分清是「问题病态」还是「算法不佳」。2. 算法设计2.1 金标准后向稳定性衡量一个数值算法有两条路线前向误差算出的答案 f̂(x) 离真值 f(x) 有多远后向误差算出的答案是不是「某个 nearby 问题的精确解」——即是否存在 x̃ 使得 f̂(x) f(x̃)且 x̃ 离 x 很近。定义算法是后向稳定的当且仅当对每个输入 x都存在 x̃ x(1 O(u)) 使得算法输出恰为 f(x̃)。直观说法它「精确地解决了一个几乎一样的问题」。为什么后向稳定能成为金标准因为它和条件数组合后给出完整的前向误差界前向相对误差≲κ(x)⋅u\text{前向相对误差} \lesssim \kappa(x) \cdot u前向相对误差≲κ(x)⋅u即后向稳定算法 × 良态问题 可靠结果。这个框架还有一个工程上极好的性质——归因清晰结果不准要么问题病态κ 大换问题建模要么算法不后向稳定换算法。为什么不直接用「前向误差小」做标准前向误差同时由问题和算法决定库作者无法对任意用户输入做出承诺而后向稳定性是算法自身的性质一次证明、处处成立。这就是工业界选择它的原因——不是它更「强」而是它可认证、可归因。2.2 灾难性抵消减法不背锅数据先脏了两个相近的浮点数相减有效数字大量丢失这就是灾难性抵消catastrophic cancellation。定量地设 x̂ x(1δ₁)、ŷ y(1δ₂)|δ| ≤ u则∣误差∣∣x^−y^∣≤u⋅∣x∣∣y∣∣x−y∣\frac{|\text{误差}|}{|\hat x - \hat y|} \le u \cdot \frac{|x| |y|}{|x - y|}∣x^−y^∣∣误差∣≤u⋅∣x−y∣∣x∣∣y∣当 x ≈ y 时放大因子 ≈ 2|x|/|x−y|可以任意大。一个关键澄清减法本身是精确的。Sterbenz 引理保证当 x/2 ≤ y ≤ 2x 时x − y 没有舍入误差。抵消之所以「灾难」是因为相减的两个数在更早的运算里已经被污染减法只是把污染暴露出来。常见场景与稳定化改写本模块实现其中三个场景抵消来源稳定化改写二次方程求根 −b±√Δb² ≫ 4ac 时 √Δ ≈ |b|q 公式 韦达定理见 2.3方差 E[x²] − (E[x])²大均值小方差时两项相近Welford 在线算法课后练习(1x) − 1x 很小1 吞掉 x 的低位log1p(x)/expm1(x)√(x1) − √xx 大时两值相近有理化1/(√(x1)√x)(f(xh)−f(x))/hh 小时分子相近复步微分/符号微分毕业项目方向 3 的伏笔二次求根的稳定公式教科书公式 x (−b ± √Δ)/(2a) 在 b 与 √Δ 同号时必然抵消。稳定做法分两步构造 q −0.5·(b sign(b)·√Δ)保证括号内同号相加、绝不相消得一个远离抵消的根 x₁ q/a另一个根用韦达定理 x₁x₂ c/a 求出x₂ c/q同样避开减法。复杂度仍为 O(1)一次 sqrt、常数次四则运算代价几乎为零误差从「灾难级」降到 O(u)——M01 实验会实测这一差距。2.3 求和算法从 O(nu) 到 O(u)n 个数朴素累加的最坏误差界约为∣误差∣≲(n−1) u∑i∣xi∣|\text{误差}| \lesssim (n-1)\, u \sum_i |x_i|∣误差∣≲(n−1)ui∑∣xi∣误差随 n线性增长且与数据顺序有关。两个改进方向Kahan / Neumaier 补偿求和。核心思想用补偿变量 c 收集每一步被舍入丢掉的部分最后加回。经典 Kahan 假设「当前和总比新加入的数大」当数据顺序打破这一假设时会失效Neumaier 改进版补上了这个角落matlib 采用的就是它。伪代码s ← 0; c ← 0 for x in data: t ← s x if |s| ≥ |x|: c ← c ((s − t) x) # x 的低位被丢进这一项 else: c ← c ((x − t) s) # s 的低位被丢进这一项 s ← t return s c复杂度 O(n)每步多常数次浮点运算误差界降到约 2u·Σ|x|与 n 基本无关。备选方案对比为什么不选它们按绝对值排序后从小到大累加直觉正确但需要 O(n log n) 排序、必须拿到全部数据不适合流式且最坏误差界仍与 n 有关。弃用但思想保留pairwise两两分组求和把序列递归二分成树形求和误差 O(log n · u)无需维护补偿变量、缓存友好、可向量化。这是NumPy 的工业选择——它是「足够好 足够快」的折中细节见第 4 节与 M12。matlib 的选择流式精确场景用 Neumaierkahan_sum大规模向量化场景留待 M09 之后再考虑 pairwise。3. 库级实现3.1 matlib 骨架本模块为 matlib 建立目录与构建骨架code/ ├── CMakeLists.txt # C17MSVC 下开启 /utf-8 /W4 ├── matlib/ │ └── float_info.hpp # M01 数值底座header-only └── test_m01.cpp # 数值验证程序22 项断言约定从本模块开始生效C17、零第三方依赖、命名向 Eigen 看齐camelCase。3.2 float_info.hpp 关键接口位级视图。bits_of用 memcpy 做标准允许的位级双关decompose按 1|11|52 拆出三段inlinestd::uint64_tbits_of(doublex){std::uint64_tb;std::memcpy(b,x,sizeof(b));// memcpy 是标准允许的位级双关写法returnb;}机器精度实测。循环对半逼近直到 1 eps/2 被舍入回 1templatetypenameTTmachine_epsilon(){T epsT(1);for(intguard0;guard256;guard){// 防御式上限if(T(1)eps/2T(1))break;eps/2;}returneps;}ULP 比较。先把位模式映射成保序整数负数区间翻转、−0 归一到 0再做无符号差——用 unsigned 是为了覆盖 ±DBL_MAX 之间距离超过 int64 范围的极端边界inlinestd::int64_tordered_bits(doublex){conststd::int64_tbstatic_caststd::int64_t(bits_of(x));return(b0)?(INT64_MIN-b):b;}inlineboolnearly_equal(doublea,doubleb,std::uint64_tmax_ulp4);边界探测。从 1.0 出发倍增直到溢出、倍减直到即将进入非规格化区实测浮点世界的两堵墙2^1023 与 2⁻¹⁰²²实现见test_m01.cpp第 [8] 段。Neumaier 补偿求和kahan_sum、稳定求根solve_quadraticq 公式 韦达定理处理 a0 退化、Δ0、二重零根三个边界、稳定 hypothypot_stable按最大值缩放。完整代码见code/matlib/float_info.hpp。3.3 数值验证实测结果测试环境Windows x64MSVC 14.44VS 2022/std:c17 /O2 /EHsc /utf-8。完整验证程序为test_m01.cpp22 项断言全部通过。摘录关键输出eps_mach(double) 2.2204460492503131e-16 # 与理论值 2^-52 一致 eps_mach(float) 1.1920929e-07 # 2^-23 bits(0.1) 0x3FB999999999999A # 无限循环二进制被截断的证据 0.10.2 0.30000000000000004 0.3 0.29999999999999999, ulp 距离 1 # 只差 1 个 ULP naive_sum 0, kahan_sum 1000 # {1e16, 1×1000, -1e16}精确和 1000 sum(0.1 x 1e6): naive |err| 1.33e-06, kahan |err| 0 x²1e8x1 的小根 参考值 -1.0000000000000002e-08 naive 求根 -7.4505805969238281e-09 相对误差 0.255 # 错了 25% 稳定版求根 -1e-08 相对误差 1.65e-16 # 回到 u 量级 sqrt(3e200²4e200²) inf溢出 hypot_stable 5e200 ((1x)-1)/x 0 log1p(x)/x 1 (x1e-16) 边界探测: 最大 2 的幂 2^1023最小正规格化数 2^-1022 2^-1022 之下还有 52 级非规格化数直到 2^-1074 才归零——gradual underflow三个值得盯着看的结果朴素求和得 0ulp(10¹⁶) 2每个 1 恰好落在舍入的中点上被丢弃——不是「误差小」是「信息整体蒸发」kahan 求和误差为 010⁶ × fl(0.1) 与 100000 的差距约 5.55×10⁻¹²小于该量级的半个 ULP约 7.3×10⁻¹²最后的加法舍入恰好把它吸收回 100000.0naive 求根错 25%教科书公式在 b 10⁸ 时小根已完全不可信——这是「算法不稳定」与「问题病态」的区分样本问题求这个根本身良态是算法把误差放了进来。4. 工业级解析4.1 EigenNumTraits 暴露机器参数Eigen 把本模块手测的所有机器参数收进一个 traits 类NumTraitsT源码坐标Eigen/src/Core/NumTraits.hEigen 3.4.0epsilon()机器精度对应 matlib 的machine_epsilondummy_precision()经验容差epsilon 的小倍数用于近似比较——对应nearly_equal的思路但 Eigen 选了「容差倍数」路线我们选了「ULP 距离」路线后者在跨量级比较时更稳highest()/lowest()取值范围边界。后续模块中 Eigen 的各类阈值判断比如秩判定都从 NumTraits 取数——机器参数集中管理是库级设计的基本功。4.2 NumPyfinfo 与 pairwise 求和NumPy 侧的对应物NumPy 2.xnp.finfo(float)返回 eps / tiny / max分别对应机器精度、最小规格化数、最大可表示数np.spacing(1.0)返回 1.0 处的 ULP恰为 2⁻⁵²np.nextafter取相邻浮点数——正是 3.2 节ulp_distance的原料np.add.reduce即np.sum的底层不用朴素累加而是pairwise 求和NumPy 1.9 起源自社区对求和精度的长期讨论不超过 128 个元素的块内做展开直加超过则递归二分。实现位于源码numpy/_core/src/umath/loops.c.srcNumPy 2.0 之前为numpy/core/src/umath/loops.c.src可在源码中 greppairwise定位M12 会逐行拆。对照 Python 生态的另外两个选择math.fsum维护全部部分和精确但慢CPython 3.12 起内置sum()改用 Neumaier 补偿——与 matlib 的选择同款。4.3 C 标准库稳定化改写早就写进了标准C99 在math.h中加入log1p、expm1、hypot、fma——这些函数存在的唯一理由就是本模块第 2 节的稳定化改写。工业共识稳定形式是默认形式朴素写法是调用者的责任。4.4 差距分析教学版实现与工业版的差距不在算法思想本模块的公式与工业界完全同源而在工程维度无 SIMD、无多精度回退、无全格式覆盖float16/bfloat16、无平台特殊值穷举测试。这些维度分别对应后续模块SIMD → M09格式与 dispatch → M10–M12。5. 随堂实验目的在自己的机器上复现本模块全部结论——机器精度不是背下来的 2⁻⁵²是测出来的。验收数据要求实测 εmachfloat/double并与理论值对照构造一组数据使朴素求和误差可见而补偿求和精确复现二次求根的 naive vs 稳定版误差差给出hypot溢出边界的实测值。结论必须回答为什么 kahan 求和在「0.1 × 10⁶」实验里误差恰好为 0提示想清楚最终那次加法的舍入。详细步骤见lab.md将随「生成 M01 实验」指令产出。6. 本模块小结double 只有 2⁶⁴ 个值误差是表示系统的固有契约不是 bugULP 随指数翻倍浮点世界必须用相对误差说话εmach 2⁻⁵²舍入单元 u 2⁻⁵³是后续所有误差分析的计量单位条件数是问题的先天体质先判良态病态再谈算法金标准是后向稳定性前向误差 ≲ κ·u归因清晰、可认证灾难性抵消中减法本身精确Sterbenz锅在上游已被污染的数据Neumaier 补偿求和把求和误差从 O(nu) 压到 O(u)NumPy 用 pairwise 走了另一条同样成立的工业路线。在主线上的位置本模块四个阶段全部走通但重心在「数学原理 算法设计」——它提供的是贯穿全课程的世界观与误差语言库级实现刚打下 matlib 的地基真正的性能与架构主题从 M08 开始。7. 课后练习概念复盘εmach 2⁻⁵² 而相对舍入误差界 u 2⁻⁵³差一个因子 2。用自己的话解释这两个数各自的定义以及为什么「测 εmach」和「用 u 做误差分析」都合理。实现给 matlib 增加dot_stable用 Neumaier 补偿的向量点积。构造两组量级悬殊的向量实测 naive 点积与dot_stable的误差差。实现实现 Welford 在线方差算法单遍、数值稳定与「E[x²] − (E[x])²」两遍公式在「均值 10⁹、标准差 1」的数据上对比复现后者算出负方差的翻车现场。挑战选做找到 NumPyloops.c.src中的 pairwise 求和实现解释块阈值为什么取 128提示想想缓存行与展开因子并复现「pairwise 误差介于 naive 与 Kahan 之间」的曲线。