浮点数不与自己一致

# 浮点数与自己不一致 来源：https://docs.merca.earth/blog/floats-dont-agree-with-themselves 我调试了一个多边形重叠测试，它在本地工作但在服务器上失败。相同的代码。相同的输入。不同的答案。决定答案的函数很小。三个点 A、B、C；返回 `(B - A) 叉积 (C - A)` 的符号。这个符号告诉你一个顶点是凸的还是凹的，一条对角线是否在多边形内部，一个点是在一条线之上还是之下。在 x86 上，LLVM 将乘法和减法合并成一个 `fma`——一次舍入而不是两次。在 WASM 上没有 FMA，所以有两次舍入。一个位于零的 epsilon 邻域内的顶点落在了相反的两侧。由于这一位的分歧，整个分解岔开了。 这不是我的代码中的错误。这是 IEEE 754 按预期工作。标准确定了存储格式。它没有确定行为。重新关联、FMA 收缩、80 位的 x87 寄存器、非规格化数刷新标志——四种不同的方式可以从相同输入得到不同答案。收紧容差也无济于事：凸分解是离散的。一个顶点要么是凹的，要么不是。输入中的 epsilon 在输出中就是二进制的差异。 所以我写了 `exact-poly`，一个完全不使用浮点数的 2D 几何库。 一个 10 顶点的星形多边形分解为 5 个凸部分，顶点标签为 v0–v9，部分标签为 p0–p5。 *一个 10 顶点的星形分解为 5 个凸部分。顶点标签（`v0`–`v9`）与输入顺序一致；部分标签（`p0.0`–`p5.3`）由级联产生。各部分面积之和等于环面积，精确相等。* 试试看： - **GitHub**—mercaearth/exact-poly (https://github.com/mercaearth/exact-poly) - **线上演示**（7 个标签页）—exact-poly.merca.earth (https://exact-poly.merca.earth/) - **生产环境**—merca.earth (https://merca.earth/) 在 WASM 客户端和验证端运行 `exact-poly`；选择一个地点，绘制一个多边形，观察相同的 `i64` 数学在两侧接受或拒绝声明 ## 为什么浮点数并非处处相同 IEEE 754 规定了格式和基本操作。这些操作之间的间隙就是可重现性泄漏的地方。 | 机制 | 效果 | 分歧之处 | |------|------|----------| | 中间寄存器 | x87 FPU 将值保留在 80 位，仅当溢出到内存时才舍入 | x86 vs ARM、RISC-V、WASM（无扩展精度） | | 融合乘加 | `fma(a, b, c)` 一次舍入而非两次——更精确，结果不同 | LLVM 在一个目标上启用，在另一个目标上禁用 | | 重新关联 | `(a + b) + c` 重写为 `a + (b + c)` 改变舍入 | `-ffast-math` 默认为开启状态，常常并非有意 | | 非规格化数 | 刷新为零避免慢速子正规数 | 每个进程的标志，被静默翻转 | “IEEE 754 兼容性”是关于浮点数如何存储，而非它们如何行为。可重现性必须明确要求：一个架构、一个工具链、一组标志。跨进程或 ISA 边界，保证就消失了。 ## 为什么一位就能破坏一个分解 `cross_sign(A, B, C)` 是算法拥有的唯一信号。正数表示左转，负数表示右转，零表示共线。从这个符号得出凹/凸、内部/外部、上方/下方。接近零时——三个点几乎共线——浮点数可能在一台机器上返回 `+1e-12` 在另一台上返回 `-2e-13`。翻转那一位符号位，不同的顶点就变成了凹的。不同的凹顶点意味着不同的起始切割，进而意味着不同的分解图。一位，多米诺骨牌效应。 “在比较之前舍入到微米”救不了你。凸分解是离散的：一个顶点要么是凹的，要么不是。任何围绕零的 epsilon 区域最终都会捕获一个真实的多边形，行为就会分叉。 诚实的修复是不使用浮点数。对于整数，符号就是位，而不是近似值。相同的输入，在 x86、ARM、WASM——任何存在 `i64` 的地方，甚至在没有浮点类型的平台上——结果相同。Shewchuk 1996 年关于自适应精度谓词的论文是 Triangle 和 CGAL 依赖的基础；每个人都在修补同一个漏洞。`exact-poly` 采取了最直接的路线：在 `i128` 中计算叉积，不进行过滤。每次调用多一个乘法器，没有“尝试双精度，失败再扩展”。始终便宜且正确。 ## 选择尺度 “不使用浮点数”很简单，但有一个昂贵的推论：整数算术不会自动缩放。你需要为你的问题选择一个尺度，并为这个选择付出代价。预算为 `i64::MAX ≈ 9.2 × 10^18`。两个旋钮：度量单位，以及你曾经想要处理的最大坐标。它们的乘积达到了上限。 | 领域 | 尺度 | 1 单位 | 最大坐标 (`i64`) | |------|------|--------|-------------------| | 测地学 (Web Mercator) | 10^6 | 1 μm | ±9 × 10^12 m — 超过行星 | | 物理引擎 | 10^3 | 1 mm | ±9 × 10^15 m — 太阳系范围 | | CAD 加工 | 10^9 | 1 nm | ±9 × 10^9 m — 工厂车间，而非城市 | 没有正确的尺度。每个都是权衡。`exact-poly` 目前为测地学校准。`SCALE = 1_000_000`。一米是一百万个单位；一个单位是一微米。消费级 GPS 在好天气下精度约为厘米，因此一微米比我们测量的任何东西都精细四个数量级——不是为了用户，而是作为边缘归一化、中点和吸附的余量。Web Mercator 中的地球周长是 `40_075_017` 米，即 `4 × 10^13` 单位。这是用例的上限，而非库的上限。从此处到 `i64::MAX` 还有五个半数量级的余量。地球上的每个坐标都适合 `i64`。 `SCALE` 和 `MAX_WORLD` 在编译时被故意硬编码。在运行时更改它们，两个运行相同代码的进程可能悄悄地对于什么是“一米”产生分歧——这正是库要防止的失败。分叉并重新编译以进行不同的校准只需要两分钟的工作；运行时配置性是自找麻烦。 ## 当 `i128` 出现时 叉积是 `(bx - ax)(cy - ay) - (by - ay)(cx - ax)`。差值最大可达 `~8 × 10^13`；它们的乘积可以达到 `6.4 × 10^27`——远远超过 `i64`。`i128::MAX` 约为 `1.7 × 10^38`，十个数量的余量。规则只有一行：坐标是 `i64`，任何相乘的都是 `i128`。 ``` // src/signed.rs:21-27 pub fn cross_sign(ax: i64, ay: i64, bx: i64, by: i64, cx: i64, cy: i64) -> i128 { let dx1 = (bx as i128) - (ax as i128); let dy1 = (by as i128) - (ay as i128); let dx2 = (cx as i128) - (ax as i128); let dy2 = (cy as i128) - (ay as i128); dx1 * dy2 - dy1 * dx2 } ``` 先扩展再减法：`(bx as i128) - (ax as i128)`，而不是 `(bx - ax) as i128`。如果先在 `i64` 中减法，乘法会溢出；你会得到 `i128` 形状的垃圾。我读过的半数天真的整数几何实现都搞反了这点。 为什么用 `i128` 而不是像 GMP 那样的任意精度有理数？因为 `i128` 适合两个寄存器，在硬件中运行，不分配任何东西。有理数需要堆和分支。对于叉积，所需的精度是有界的：输入 `i64`，两个差值的乘积，适合 `i128`。不是无限，只是足以满足一个经过证明的界限。与 Shewchuk 在过滤谓词中的技巧相同：便宜且正确优先，仅在必须时才回退。[1] ## 面积作为精确等式 鞋带公式给出 `2 × 面积`。没有理由除以 2——每个检查都在加倍后的值上进行，而除法会丢弃低位。为什么要丢弃一个你用 `i128` 付出的位？`twice_area` 是一个有符号的 `i128`。正数表示逆时针，负数表示顺时针。库的核心不变式是： ``` sum(twice_area(parts[i])) == twice_area(original_ring) ``` **`==`，而不是 `abs(a - b) < eps`。** 精确整数等式。对于浮点数，这个不变式在原则上是不可能的——累加 N 个面积会以不同方式扰动低位。对于整数，各部分加起来等于整体，逐位相等。如果任何一位不对，分解器就丢失或复制了一个微观碎片，测试就会失败。 ## 为什么还要分解 分离轴定理是一个快速、简单的测试，用于判断两个多边形是否重叠。它只适用于凸多边形。任何非凸多边形都必须被切割成凸部分然后成对检查。最简单的非凸形状是 L 形。一个凹顶点就让 SAT 无法适用。一个典型的地域多边形——一块地块、一个频段、一个传感器覆盖区域——有三到五十个顶点和几个凹角。SAT 需要凸部分才能诚实回答“这些是否重叠”。 ## 带旋转的级联 凸分解的算法自八十年代以来就存在。Chazelle 在 1991 年发表了渐近最优的 `O(n + r^2)`，这个上限至今仍然成立。用整数算术诚实地实现它需要三个月的工作和一个不稳定的梯形分解来保持直立。对于实际的 `n < 100`，更简单的算法比 Chazelle 的实现运行得更快。教科书中有 Chazelle 1991。生产环境中有 Bayazit、Keil 和 Hertel-Mehlhorn。 每个更简单的算法都有弱点。ExactPartition 是一个贪心分区器，优先考虑来自凹顶点的对角线——快速，通常有效，但在病态输入上会放弃。Bayazit 是递归的，可以在边的中点插入 Steiner 点——更强，但未必能在合理的步数内收敛。EarClip 是 1985 年的耳剪三角剖分；它在任何简单多边形上工作并产生 `O(n)` 个三角形，这些三角形随后需要合并。这就是 Hertel-Mehlhorn (1983) 登场的地方——一个合并步骤，提供了 `4 × OPT` 的部分数量保证。 诀窍不是选择一个。而是将它们串联起来。如果第一个失败，下一个接手。如果所有三个都失败，旋转环再试一次。纵深防御。 *图 1. 级联。每个算法尝试；失败时下一个接手。如果所有三个都失败，旋转环并重试。* Hertel-Mehlhorn 合并步骤本身就是整数如何使事情变得便宜的一个小例子。要合并两个共享边的三角形，你需要检查结果是否保持凸性：两个符号相同的叉积。没有 epsilon，没有容差。 旋转步骤听起来像是“再试一次并祈祷”。这是一个启发式——没有最优性证明。它之所以有效，是因为病态通常是一个顶点的函数，而不是所有顶点。起始索引是一个你可以拉动的杠杆，而不改变几何形状：相同的面积、相同的周长、相同的凹顶点集合。不同的遍历，不同的结果。最坏情况是 `N^2`；实际上，几乎一切都在旋转零时解决。 ## 一个接缝处的三个陷阱 Bayazit 有时会在边的中点插入一个顶点——一个 Steiner 点——当没有对角线能在现有顶点之间工作时。很有用，但它点亮了三个整数特定的陷阱。 ### 不丢失位的中点 `(a + b) / 2` 在整数上会丢掉低位；点最终会偏离半个单位，与邻居的共线性会破坏。修复是 `round_div2`——`(a + b) >> 1`，对于负数使用远离零的舍入。 ### 远离端点吸附 如果中点靠近一个现有顶点，它必须被推开，否则三个点会共线并破坏通过它们的每条对角线的 `in_cone`。 ### 合成顶点记账 Steiner 点会被标记，这样它们不会改变输入顶点计数。它们在分解期间内部存在，仅在提交时具体化为外环。否则，下游验证器会计数到与客户端创建的顶点数不同的数量，系统再次默默地与自己不一致——这正是整数重写本应消除的失败模式。 这些不是天才的想法。它们是三个小楔子，打入同一个接缝——连续几何与离散算术相遇之处。文献中大多跳过它们。当针对特定多边形的测试拒绝收敛时，它们就会出现。 ## 什么不起作用 诚实的清单，不是营销清单。 - **自相交输入**——被 `is_simple()` 拒绝（`O(n^2)` 检查非相邻边对）。一个 8 字形返回 `DecompError::NotSimple`。没有英雄式恢复；在调用者端修复它，例如使用 Clipper2 布尔预处理。 - **带孔的多边形**——不支持。仅外环。 - **部分超过 `MAX_VERTICES_PER_PART`**——验证器拒绝。不是回退；是信号表明输入不适合配置的限制。 - **零面积或共线三顶点多边形**——在 `normalize_ring` 折叠连续共线顶点后被拒绝。 - **`n > 1000` 的基准测试**——未发布。对于 `n < 100`，级联运行在毫秒级（主观——我没有正确运行 `criterion`）。当 `n = 1000+` 时，旋转 × `is_simple` × 算法成本可能达到秒级。 库不试图解决的问题，它在输入时拒绝。没有未定义行为，没有关于 8 字形的无限循环。 ## 回到确定性 级联是一个条件分支树。每个分支归结为叉积的符号、`i128` 面积的比较或 `i64` 长度的比较。任何地方都没有浮点数。相同的环在所有架构上产生相同的分解：在浏览器中（`wasm32`）、在服务器上（`x86_64`）以及在根本没有浮点类型的运行时——一个没有 FPU 的 Cortex-M0、一个 ZK 电路 IR、一个 Move VM。一个 Rust crate，多个编译目标，处处逐位相同的结果。[2] *图 2. 一个源 crate，五个编译目标，处处相同的 `i64`/`i128` 输出。* ## 测试作为不变式 存储库的测试围绕属性组织，而非覆盖率。全局的是面积守恒：在任何成功分解后，部分面积之和等于原始环面积，精确相等。一个断言就能捕获几乎任何局部错误。 ``` let original = twice_area_fp2(&ring); let result = decompose(&ring, algo, cfg)?; let sum: i128 = result.parts.iter() .map(|p| twice_area_fp2(p)) .sum(); assert_eq!(sum, original); ``` 回归测试携带它们的历史： | 测试属性 | 保护内容 | |----------|----------| | `l_shape` | 典型凹案例——如果它坏了，一切都坏了 | | `comb_twelve_vertices` | 十二顶点梳形——早期版本分类错误 | | `sat_touching_edges_do_not_overlap` | 严格 `>` 重叠，而非 `>=`——共享边的正方形是接触而非重叠 | | `cross_sign_no_overflow_with_max_world_coordinates` | 叉积在 `4 × 10^13` 处落于 `1.6 × 10^27`，完全在 `i128` 内 | 最后一个存在是为了让下一个想知道“我们能降级到 `i64` 吗？”的人看到一个失败的测试。 `validate_onchain.rs` 镜像了验证器的检查：部分数量、顶点数量、最小边长、紧密度限制。相同的逻辑，相同的步骤顺序。如果 `decompose` 和 `validate_onchain` 在本地通过，它们在远程也会通过。不是“可能”。总是——代码在两边是同一个二进制。 ## 结语 浮点数关乎精度。整数算术关乎一致性。当两个进程运行相同的算法并且必须得到相同的答案时，浮点数是敌人。`exact-poly` 放弃了精度以获得确定性。不是更精确，而是完全精确。不是更快，而是在任何有 128 位整数的地方都相同。如果你想要 IEEE 754 兼容性，有几十个库可用。如果你想要逐位相同的结果，只有一个。 [1] Shewchuk, Jonathan Richard. “Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates.” Discrete & Computational Geometry 18, no. 3 (1997): 305–63. [2] 在撰写本文时。我还没有在 ZK 电路上运行它。这是合理的目标架构列表，而非性能声明。