几何谓词orient3d讲解

引言：浮点数与几何计算的冲突

在计算机图形学、计算几何、CAD建模和物理模拟等领域，我们不断地需要程序做出基础的几何判断。例如：

• 三个点 A, B, C 是顺时针排列还是逆时针排列？(orient2d)
• 一个点 D 是在一个由 A, B, C 构成的三角形内部还是外部？
• 一个点 D 是在由 A, B, C 定义的平面之上还是之下？(orient3d)
• 一个点 E 是在由 A, B, C, D 定义的球体内部还是外部？(insphere)

这些判断被称为几何谓词 (Geometric Predicates)。它们是构建几乎所有复杂几何算法（如构建三角网、计算凸包、布尔运算）的基石。在理论上，这些判断是通过计算行列式的符号来完成的。例如，orient3d 的正负号取决于下面这个行列式：

$$\text{orient3d}(A,B,C,D)=\text{sign}\det \left(\begin{array}{ccc}{A}_x-D_x& A_y-D_y& A_z-D_z\\ B_x-D_x& B_y-D_y& B_z-D_z\\ C_x-D_x& C_y-D_y& C_z-D_z\end{array}\right)$$

然而，当理论撞上计算机硬件的现实时，一场“看不见的危机”便发生了。我们几乎所有的计算都依赖于硬件实现的浮点数（如64位的 double）。浮点数虽然速度极快，但其精度是有限的。它无法精确表示所有实数，计算过程中会产生微小的舍入误差 (Rounding Error)。

对于大多数计算，这点误差无关紧要。但在几何谓词中，它却是致命的。当输入的点共线、共面或共球（或极其接近这些“退化”状态）时，行列式的真实值会非常接近于零。此时，一个微不足道的舍入误差就可能导致计算结果的符号从正变为负，或从非零变为零。这种错误的符号判断会像病毒一样在算法中传播，导致程序产生荒谬的结果，比如生成一个扭曲的网格、错误的拓扑结构，甚至陷入无限循环而崩溃。

核心问题：如何设计一个既能利用硬件浮点数的超高速度，又能保证在所有情况下（尤其是退化情况下）都做出100%正确几何判断的谓词？

方法论：Jonathan Shewchuk 的自适应精度框架

1997年，Jonathan Richard Shewchuk 在其里程碑式的论文中提出了一个极其优雅且高效的解决方案，至今仍是该领域的黄金标准。其核心思想是：为绝大多数“简单”情况提供极速的答案，并为极少数“困难”情况提供绝对正确的答案。

这个框架被称为自适应精度浮点算术 (Adaptive Precision Floating-Point Arithmetic)，它通过一个多级过滤系统来实现：

1. 静态过滤器 (Static Filter): 这是第一道防线。

   • 计算：首先，使用标准的、速度飞快的 double 类型直接计算行列式的近似值 det_approx。
   • 误差分析：然后，通过一个基于输入的坐标值、预先计算好的、绝对可靠的公式，估算出这次计算可能产生的最大舍入误差 error_bound。这个误差界是经过严格数学证明的，保证涵盖了最坏情况。
   • 判断：比较 det_approx 和 error_bound。如果 |det_approx| > error_bound，这意味着近似值的绝对值远大于可能产生的最大误差。因此，舍入误差绝不可能改变结果的符号。此时，我们可以确信 det_approx 的符号就是正确的符号，直接返回即可。

   在实践中，超过99.9%的计算都会在这里通过，我们以接近硬件的速度得到了可靠的结果。

2. 精确算术 (Exact Arithmetic): 当过滤器失败时（即 |det_approx| <= error_bound），程序便知道不能相信这个结果。此时，它会无缝切换到第二道防线：一个用纯软件模拟的、无误差的高精度算术系统。这个系统本身也是“自适应”的，它只引入必要的精度来解决问题，从而将性能开销降到最低。

数学核心：无误差变换与浮点数展开式

这个框架的魔力源于一个深刻的数学事实：两个 double 浮点数进行一次基础算术运算（加、减、乘）的真实数学结果，可以被无损地表示为另外两个 double 的和。

1. 原子操作 (Error-Free Transformations)

Shewchuk 的代码实现了一系列底层的“原子”操作，它们是构建精确计算大厦的砖块：

• Two-Sum(a, b) -> (x, y): 这个操作计算 x = a + b（硬件浮点加法），并精确地捕获其舍入误差 y = (a + b) - x。最终，a + b 在数学上严格等于 x + y。
• Two-Prod(a, b) -> (p, e): 类似地，它计算 p = a * b，并精确捕获其舍入误差 e = (a * b) - p。这背后的原理是：一个53位有效数字的 double 乘以另一个53位的 double，其精确结果最多需要106位来表示。这106位的结果可以完美地被拆分成高位的53位（存入 p）和低位的53位（存入 e）。

这些操作被称为无误差变换，因为它们没有丢失任何信息，只是将一次运算的结果信息分布到了两个 double 中。

2. 数据结构：浮点数展开式 (Floating-Point Expansion)

有了原子操作，我们就可以定义一个高精度数字的表示方法：浮点数展开式。
一个展开式就是一个 double 数组，例如 E = [d1, d2, d3, ..., dn]。它所代表的数学值就是 d1 + d2 + d3 + ... + dn。

这个数组有一个至关重要的性质：非重叠 (Non-overlapping)。这意味着每个分量 di 的值都远小于它前一个分量 d(i-1) 的最低有效位。这保证了在简单相加时，它们的尾数不会重叠，从而不会引入新的舍入误差。

基于这个数据结构，可以定义一整套高精度算术库：

• expansion_sum(E, F): 将两个展开式相加，生成一个更长的展开式。
• scale_expansion(E, b): 将一个展开式乘以一个 double，生成一个新的展开式。

3. orient3d 的精确计算路径

当 orient3d 的快速过滤器失败后，它的精确计算路径 (orient3dadapt) 就会启动：

1. 第一轮精确计算：它不再将坐标差 adx, bdx 等看作是 double，而是将它们看作是长度为1的展开式。行列式中的每一个乘法都用 Two-Prod 代替，每一次加减法都用 expansion_sum 代替。这会产生一个代表行列式精确值的展开式 D_exact。

2. 自适应的本质：这个过程不是一步到位的。算法会先使用较短的展开式计算出一个中间结果，然后再次进行过滤判断。如果这个更高精度的结果已经足以确定符号，就立即返回。只有在符号仍然不确定的极端情况下，它才会引入更多的误差项（例如，初始坐标减法 pa[0] - pd[0] 本身的误差），将展开式变得更长，直到最终的符号水落石出。

结果与影响

当整个精确计算完成后，我们会得到最终的展开式 D_exact = [d1, d2, ..., dn]。如何得到最终的符号？

• 不是将它们再次相加（这会重新引入误差）。
• 而是直接返回展开式中最重要分量的符号。由于 Shewchuk 的展开式是按绝对值排序的，这通常是最后一个非零元素 dn 的符号。这个符号保证是数学上完全正确的。

总结该框架的成果：

1. 健壮性 (Robustness): 这是最大的成就。采用此方法的几何算法永远不会因为浮点数舍入误差而失败。它们可以正确处理任何退化或近退化的输入，这对于需要高可靠性的工业软件（CAD/CAM, GIS）至关重要。

2. 高性能 (Performance): 由于绝大多数情况都被极快的硬件过滤器处理，其平均性能非常接近于完全不考虑健壮性的“天真”算法。只有在极少数情况下，才需要付出精确计算的代价。我们得到了两全其美的结果：平时的速度和关键时刻的正确性。

3. 行业标准: Shewchuk 的论文和其公开的 predicates.c 代码已经成为计算几何领域的基石和事实上的标准。高质量的计算几何库（如 CGAL, Geogram）的核心都构建在这一理论之上。

4. 现代发展: 尽管核心思想保持不变，但现代硬件的发展进一步提升了其性能。例如，现代CPU中的 FMA (Fused Multiply-Add) 指令可以一步完成 a*b+c 并且只有一次舍入，这使得 Two-Prod 的实现变得极其简单和高效，从而大幅加速了精确计算路径。对于大规模问题，这种过滤方法也非常适合在 GPU 上进行并行化处理，实现数量级的性能飞跃。