1. 浮点数为什么反直觉
在实数域里 (a + b) + c == a + (b + c) 天经地义，但在浮点世界里它可能不成立。原因很简单：浮点数只有有限位，每一次运算都要舍入。理解数值计算的第一课，就是放弃「浮点等于实数」的直觉，转而用误差的眼光看待每一步操作。
>>> 0.1 + 0.2 == 0.3
False
>>> 0.1 + 0.2
0.30000000000000004
>>> 1e16 + 1 == 1e16
True # 1 被"吞掉"了
这三行输出分别对应：表示误差（0.1 无法精确表示）、舍入误差累积、大数吃小数。它们不是 bug，而是有限精度的必然结果。
2. IEEE 754 表示
2.1 位布局
IEEE 754 把浮点数编码为「符号位 + 指数位 + 尾数位」：
| 格式 | 总位宽 | 符号 | 指数 | 尾数 | 偏置 |
|---|---|---|---|---|---|
| binary32（单精度） | 32 | 1 | 8 | 23 | 127 |
| binary64（双精度） | 64 | 1 | 11 | 52 | 1023 |
数值公式（规格化数）：
value = (-1)^sign × 1.fraction × 2^(exponent - bias)
注意尾数隐含了一个前导 1（implicit leading bit），因此双精度有效位实际是 53 位（52 位显式 + 1 位隐含）。
2.2 特殊值编码
指数位全 0 或全 1 时是特殊值：
| 指数 | 尾数 | 含义 |
|---|---|---|
| 全 0 | 0 | ±0 |
| 全 0 | 非 0 | 次正规数（subnormal），无隐含 1 |
| 全 1 | 0 | ±∞ |
| 全 1 | 非 0 | NaN（quiet / signaling） |
次正规数的存在是为了渐进下溢：最小的规格化双精度数约 2.2e-308，次正规数把可表示范围延伸到约 4.9e-324，代价是精度随指数下降而逐渐损失。
// 用位运算查看浮点的原始编码
#include <stdio.h>
#include <string.h>
#include <stdint.h>
void dump(float f) {
uint32_t bits;
memcpy(&bits, &f, sizeof bits);
printf("%-12g 0x%08X sign=%u exp=%u frac=0x%06X\n",
f, bits, bits >> 31, (bits >> 23) & 0xFF, bits & 0x7FFFFF);
}
int main(void) {
dump(1.0f); // 0x3F800000 exp=127 frac=0
dump(-2.5f); // 0xC0200000
dump(0.0f); // 0x00000000
dump(1.0f/0.0f); // 0x7F800000 +Inf
dump(0.0f/0.0f); // 0x7FC00000 NaN
return 0;
}
3. 误差、机器精度与 ULP
3.1 三种误差来源
- 表示误差（Representation Error）：十进制小数转二进制时的截断，如 0.1。
- 舍入误差（Rounding Error）：每次运算结果落在两个可表示数之间，需按舍入规则取整。
- 截断误差（Truncation Error）：算法本身的近似，如用有限项泰勒级数逼近函数。
3.2 机器精度与 ULP
机器精度（machine epsilon） 定义为 1.0 与下一个可表示数之差，即 2^-52 ≈ 2.22e-16（双精度）。
ULP（Unit in the Last Place） 是「相邻两个浮点数之间的距离」，它随数值量级变化：
import math, sys
eps = sys.float_info.epsilon
print(eps) # 2.220446049250313e-16
def ulp(x):
# 返回 x 处的 ULP
if x == 0: return 5e-324
return math.ulp(x)
for v in (1.0, 1e6, 1e16, 1e-8):
print(f"x={v:<10g} ulp={ulp(v):.3e}")
# x=1 ulp=2.220e-16
# x=1e6 ulp=1.164e-10
# x=1e16 ulp=2.000e+00
# x=1e-8 ulp=1.694e-24
关键认知：误差是相对的，不是绝对的。
1e16处的 ULP 已经是 2，说明在该量级上+1根本无意义。
3.3 舍入模式
IEEE 754 定义了五种舍入模式，默认是就近舍入、平局取偶（round-to-nearest, ties-to-even）：
| 模式 | 说明 |
|---|---|
| round-to-nearest-even | 默认，误差最小 |
| round-toward-zero | 截断 |
| round-toward-+∞ | 向上 |
| round-toward–∞ | 向下 |
| round-to-nearest-away | 平局远离零 |
默认模式把「平局」时的尾数最低位取偶，长期来看能抵消系统性偏差。
4. 灾难性抵消
灾难性抵消（Catastrophic Cancellation） 发生在两个相近的大数相减时：结果的有效位数急剧减少，原有误差被放大。
4.1 经典例子
# 双精度下 1e16 + 1 - 1e16 应为 1，实际为 0
>>> (1e16 + 1.0) - 1e16
0.0
# 用直接公式解 x² - 1e8·x + 1 = 0 的两个根
>>> import math
>>> a, b, c = 1.0, -1e8, 1.0
>>> d = math.sqrt(b*b - 4*a*c)
>>> x1 = (-b + d) / (2*a)
>>> x2 = (-b - d) / (2*a)
>>> x1, x2
(1e8, 0.0) # x2 本应是 1e-8，被抵消成 0
x2 的真实值是 1e-8，但 -b - d 中两个几乎相等的大数相减，有效位损失殆尽。
4.2 稳定化改造
利用韦达定理 x1 · x2 = c / a，改用乘法求小根：
def solve_quadratic(a, b, c):
d = math.sqrt(b*b - 4*a*c)
# 先算不会抵消的那个根
if b >= 0:
x1 = (-b - d) / (2*a)
else:
x1 = (-b + d) / (2*a)
x2 = c / (a * x1) # 韦达定理，避免相减
return x1, x2
print(solve_quadratic(1.0, -1e8, 1.0)) # (1e8, 1e-8)
4.3 三种常见抵消场景
| 场景 | 表现 | 对策 |
|---|---|---|
| 相近大数相减 | 有效位丢失 | 代数变形、韦达定理 |
| 级数交错求和 | 部分和震荡 | Kahan 求和、重排 |
| 小量加到大数 | 小量被吞 | 先排序再累加、分块求和 |
5. 数值稳定性
一个算法若输入的微小扰动只引起输出的微小变化，就是数值稳定的；反之不稳定。注意稳定性与精度是两回事：稳定算法可以有误差，但不稳定算法会放大误差。
5.1 求和：朴素 vs Kahan
朴素累加在「大量小量加到一个大和」时会持续丢精度。Kahan 求和用一个补偿变量记录丢失的低位：
def naive_sum(xs):
s = 0.0
for x in xs:
s += x
return s
def kahan_sum(xs):
s = 0.0
c = 0.0 # 补偿项
for x in xs:
y = x - c # 先减去上次丢失的部分
t = s + y
c = (t - s) - y # 本次丢失的低位
s = t
return s
xs = [1e16, 1.0, -1e16] * 1000
print(naive_sum(xs)) # 0.0（每次 1.0 都被吞）
print(kahan_sum(xs)) # 1000.0
5.2 递推的稳定性
某些递推公式对初值极其敏感。例如计算积分 I_n = ∫₀¹ xⁿ/(x+5) dx，其精确递推为 I_n = 1/n - 5·I_{n-1}：
# 正向递推：误差每步放大 5 倍，很快发散
I = 0.1823215567939546 # I0 = ln(6/5)
for n in range(1, 20):
I = 1/n - 5*I
print(I) # n=19 时已严重偏离真值
# 反向递推：从大 n 向小 n，误差每步缩小 5 倍，稳定
I = 0.0
for n in range(60, 0, -1):
I = (1/n - I) / 5
print(I) # n=1 时 ≈ 0.0883922...
递推的稳定性由误差传递因子决定：
|因子| > 1则误差放大（不稳定），< 1则衰减（稳定）。方向选错，再高的精度也救不回来。
5.3 避免「相等」比较
浮点数不应直接用 == 比较，而应判断差值是否小于容差：
import math
def almost_equal(a, b, rel_tol=1e-9, abs_tol=1e-12):
return math.isclose(a, b, rel_tol=rel_tol, abs_tol=abs_tol)
# 反例：用绝对容差在大量级上失效
def bad_eq(a, b):
return abs(a - b) < 1e-9 # 1e16 附近完全失效
6. 条件数与病态问题
6.1 条件数定义
问题的**条件数（condition number）**衡量「输入相对误差」到「输出相对误差」的放大倍数：
κ = |相对输出变化| / |相对输入变化|
κ ≈ 1：良态（well-conditioned），输入小扰动 → 输出小变化。κ ≫ 1：病态（ill-conditioned），微小输入误差被剧烈放大。
# 解线性方程组 Ax = b 的条件数
import numpy as np
A = np.array([[1.0, 1.0],
[1.0, 1.0000000001]]) # 两行几乎相同
print(np.linalg.cond(A)) # 约 4e10，病态
6.2 条件数 vs 稳定性
| 概念 | 归属 | 含义 |
|---|---|---|
| 条件数 | 问题的属性 | 问题本身对扰动的敏感度 |
| 稳定性 | 算法的属性 | 算法是否放大误差 |
关键结论：病态问题无法用稳定算法拯救——误差下界由条件数决定；而良态问题若用不稳定算法，照样会算错。因此先判断问题条件数，再选择算法。
7. 迭代法与收敛
7.1 二分法
二分法（Bisection）利用介值定理，每步把区间减半，线性收敛且必定收敛：
def bisect(f, lo, hi, tol=1e-12, max_iter=200):
flo = f(lo)
for _ in range(max_iter):
mid = 0.5 * (lo + hi)
fmid = f(mid)
if fmid == 0 or (hi - lo) / 2 < tol:
return mid
if (flo < 0) != (fmid < 0):
hi = mid
else:
lo, flo = mid, fmid
return 0.5 * (lo + hi)
# 求 x³ - 2 = 0，即 ∛2
print(bisect(lambda x: x**3 - 2, 0, 2)) # ≈ 1.2599210498948732
每次迭代误差减半：e_{n+1} ≈ 0.5 · e_n，约需 log₂((b-a)/tol) 步。
7.2 牛顿法
牛顿法（Newton-Raphson）用切线逼近根，二次收敛（误差平方级下降），但需导数且可能发散：
x_{n+1} = x_n - f(x_n) / f'(x_n)
def newton(f, df, x0, tol=1e-12, max_iter=100):
x = x0
for _ in range(max_iter):
fx = f(x)
if abs(fx) < tol:
break
x = x - fx / df(x)
return x
# ∛2：f = x³-2, f' = 3x²
print(newton(lambda x: x**3-2, lambda x: 3*x*x, 1.0)) # ≈ 1.2599210498948732
7.3 收敛阶对比
| 方法 | 收敛阶 | 每步函数求值 | 优点 | 缺点 |
|---|---|---|---|---|
| 二分法 | 1（线性） | 1 | 必定收敛 | 慢 |
| 牛顿法 | 2（二次） | 2（含导数） | 极快 | 需导数、可能发散 |
| 割线法 | ≈1.618 | 1 | 免导数 | 略慢于牛顿 |
| 不动点迭代 | 视 g’ 而定 | 1 | 简单 | 收敛条件苛刻 |
7.4 终止准则
迭代不能只看 |f(x)| < tol：当根处斜率极大时，f(x) 已很小但 x 仍远离真值。稳健的终止条件应同时满足：
converged = (abs(x - x_prev) < xtol) or (abs(f(x)) < ftol)
并设置 max_iter 兜底，防止死循环。
8. 实践建议
- 能用高精度就用高精度：科学计算默认
float64；金融金额用定点数或decimal。 - 警惕相减：看到两个相近量相减，先想代数变形。
- 求和先排序：从小到大累加能显著降低舍入误差。
- 用
math.fsum/ Kahan：Python 的math.fsum内部做了精确求和。 - 比较用
isclose：不要用==。 - 检查条件数：矩阵求解前先看
cond，必要时用正则化或更高精度。
import math
# fsum 提供精确的浮点求和
print(math.fsum([0.1]*10)) # 1.0
9. 小结
浮点误差的根源是有限位表示，其表现集中在三处：表示误差（0.1 无法精确表示）、舍入误差（每次运算都要取整）、以及灾难性抵消（相近大数相减导致有效位骤减）。数值稳定性是算法的属性，条件数是问题的属性，两者共同决定最终精度。工程上，优先做代数变形、用 Kahan 求和、以 isclose 代替相等判断，并在迭代中同时约束自变量与函数值的变化。浮点编码的位级细节最终要落回计算机组成的二进制表示与浮点运算单元。
参考文章
- 计算机组成与数据表示：https://plumephp.com/cs-computer-organization/
- 位运算与整数表示：https://plumephp.com/cs-bitwise-math/
- 算法复杂度与收敛分析：https://plumephp.com/cs-algorithm-complexity/
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。