健壮的线性插值

Question

给定二维空间中两个线段端点 A 和 B，我希望能够基于一个值 t 进行线性插值，即：

C = A + t(B-A)

在理想的情况下，A、B和C应该共线。然而，由于我们使用有限的浮点数，所以会出现一些小偏差。为了解决与其他操作相关的数字问题，我使用了Jonathan Shewchuk最初创建的强大的自适应例程。具体来说，Shewchuk实现了一个方向函数orient2d，它使用自适应精度来准确测试三个点的方向。

我的问题是：是否有已知的程序可以使用浮点数计算插值，使其恰好位于A和B之间的直线上？在这里，我更关心插值本身的精度而不是结果的共线性。换句话说，只要满足共线性，C稍微偏移一下也没关系。

- MrMobster

---

如果共线性比准确性更重要，让C等于A。否则，放弃这个想法。 - user1196549

为什么不计算 C=A+t.(B-A)，然后在C周围搜索区域，选择最佳的C，使得dot(C-A,B-A)/(|C-A|.|B-A|)最接近于1。您还可以尝试cross(C-A,B-A)最小（三角形面积最小）。对于这个计算，您可以使用每个值的2个双精度浮点数来增强精度，而无需将其用于所有点... - Spektre

@Spektre，那是我的想法，使用nextafter（）在C周围搜索，直到orient（A，C'，B）== 0。但它可能会非常慢，而且我确信还有很多意外的边缘情况。恐怕Yves的帖子可能是正确的。我只是想知道是否有一些论文处理这个问题（因为我读过的大多数都从一开始就不考虑这样的问题）。 - MrMobster

@MrMobster 我认为速度不会太慢。您无需搜索大区域，只需在C周围的圆形/正方形中搜索几个坐标的ulp大小即可。为了提高精度，您还可以使用相对坐标，使点(0,0,0)成为A,B或(A+B)/2，您会惊讶于它能做到多少，请参见光线和椭球体交点精度改进。我根本看不到任何边缘情况，但是是的，您不能期望完全匹配，只能得到最佳匹配。 - Spektre

@MrMobster 顺便说一下，你可以尝试在固定点“浮点数”上执行修改后的整数 DDA，以实现“ulp”精度（不改变它们的指数）。它只使用 +,-,<，但这将需要 for 循环，因此速度会较慢，为 O(|AC|)。 - Spektre

Answer 1

坏消息

请求无法满足。对于某些A和B的值，除了0和1之外，没有其他值t使得lerp(A,B,t)为浮点数。

在单精度下的一个简单例子是x1 = 12345678.f和x2 = 12345679.f。无论y1和y2的值如何，所需结果必须具有介于12345678.f和12345679.f之间的x分量，而这两者之间没有单精度浮点数。

（有点）好消息

然而，确切的插值值可以表示为5个浮点值（在2D情况下为向量）的总和：一个用于公式的结果，一个用于每个操作的误差[1]，以及一个将误差乘以t的值。我不确定这对您是否有用。以下是单精度中使用融合乘加计算产品误差的算法的1D C版本，以便简化：

#include <math.h>

float exact_sum(float a, float b, float *err)
{
    float sum = a + b;
    float z = sum - a;
    *err = a - (sum - z) + (b - z);
    return sum;
}

float exact_mul(float a, float b, float *err)
{
    float prod = a * b;
    *err = fmaf(a, b, -prod);
    return prod;
}

float exact_lerp(float A, float B, float t,
                 float *err1, float *err2, float *err3, float *err4)
{
    float diff = exact_sum(B, -A, err1);
    float prod = exact_mul(diff, t, err2);
    *err1 = exact_mul(*err1, t, err4);
    return exact_sum(A, prod, err3);
}

为了使该算法正常工作，操作需要符合IEEE-754的最近舍入模式。这在C标准中并不保证，但GNU gcc编译器可以被指示在支持SSE2的处理器上执行此操作[2][3]。

可以保证(result + err1 + err2 + err3 + err4)的算术加法等于所需的结果；然而，不能保证这些量的浮点加法是精确的。

使用上面的例子，exact_lerp(12345678.f, 12345679.f, 0.300000011920928955078125f, &err1, &err2, &err3, &err4)返回结果12345678.f，而err1、err2、err3和err4分别为0.0f、0.0f、0.300000011920928955078125f和0.0f。确实，正确的结果是12345678.300000011920928955078125，它无法用单精度浮点数表示。

更复杂的例子：exact_lerp(0.23456789553165435791015625f, 7.345678806304931640625f, 0.300000011920928955078125f, &err1, &err2, &err3, &err4)返回2.3679010868072509765625f，错误分别为6.7055225372314453125e-08f、8.4771045294473879039287567138671875e-08f、1.490116119384765625e-08f和2.66453525910037569701671600341796875e-15f。这些数字加起来等于精确结果，即2.36790125353468550173374751466326415538787841796875，不能在单精度浮点数中精确存储。

上述示例中的所有数字都使用其精确值而非近似值编写。例如，0.3无法用单精度浮点数精确表示；最接近的一个具有精确值0.300000011920928955078125，这是我使用的那个。

也许可以尝试计算err1 + err2 + err3 + err4 + result（按照此顺序），以获得在您的用例中被认为共线的近似值。

参考资料

• [1] Graillat, Stef (2007). 《精确浮点乘法和指数运算》
• [2] 启用严格浮点模式(GCC)
• [3] GCC中浮点数运算的语义

- Pedro Gimeno