稳健的线性插值
Robust linear interpolation
给定两个线段端点A和B(二维),我想根据值t进行线性插值,即:
C = A + t(B-A)
在理想世界中,A、B 和 C 应该是共线的。然而,我们在这里使用有限的浮点运算,所以会有小的偏差。为了解决其他操作的数值问题,我使用了最初由 Jonathan Shewchuk 创建的强大的自适应例程。特别是,Shewchuk 实现了一个方向函数 orient2d
,它使用自适应精度来精确测试三个点的方向。
这是我的问题:是否有一个已知的程序可以如何使用浮点数学计算插值,以便它正好位于 A 和 B 之间的线上?在这里,我不太关心插值本身的准确性,而更关心由此产生的共线性。换句话说,只要满足共线性,C 稍微移动一点就可以了。
坏消息
无法满足请求。有 A
和 B
的值,除了 0 和 1 之外没有 t
的值,lerp(A, B, t)
是浮点数。
单精度的一个简单示例是 x1 = 12345678.f
和 x2 = 12345679.f
。无论y1
和y2
的取值如何,要求的结果必须在12345678.f
和12345679.f
之间有一个x
分量,并且在12345678.f
和12345679.f
之间没有单精度浮点数这两个。
(有点)好消息
然而,精确的插值可以表示为 5 个浮点值(在二维情况下为向量)的总和:一个用于公式的结果,一个用于每个操作中的错误 [1] 和一个用于将误差乘以 t
。我不确定这是否对您有用。下面是单精度算法的一维 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)
returns 12345678.f
和 err1
、err2
、err3
和 [= 的结果29=]分别是0.0f
、0.0f
、0.300000011920928955078125f
和0.0f
。事实上,正确的结果是 12345678.300000011920928955078125 不能表示为单精度浮点数。
一个更复杂的例子:exact_lerp(0.23456789553165435791015625f, 7.345678806304931640625f, 0.300000011920928955078125f, &err1, &err2, &err3, &err4)
returns 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] 格雷拉特,斯特夫 (2007)。 Accurate Floating Point Product and Exponentiation.
- [2] Enabling strict floating point mode in GCC
- [3]Semantics of Floating Point Math in GCC
给定两个线段端点A和B(二维),我想根据值t进行线性插值,即:
C = A + t(B-A)
在理想世界中,A、B 和 C 应该是共线的。然而,我们在这里使用有限的浮点运算,所以会有小的偏差。为了解决其他操作的数值问题,我使用了最初由 Jonathan Shewchuk 创建的强大的自适应例程。特别是,Shewchuk 实现了一个方向函数 orient2d
,它使用自适应精度来精确测试三个点的方向。
这是我的问题:是否有一个已知的程序可以如何使用浮点数学计算插值,以便它正好位于 A 和 B 之间的线上?在这里,我不太关心插值本身的准确性,而更关心由此产生的共线性。换句话说,只要满足共线性,C 稍微移动一点就可以了。
坏消息
无法满足请求。有 A
和 B
的值,除了 0 和 1 之外没有 t
的值,lerp(A, B, t)
是浮点数。
单精度的一个简单示例是 x1 = 12345678.f
和 x2 = 12345679.f
。无论y1
和y2
的取值如何,要求的结果必须在12345678.f
和12345679.f
之间有一个x
分量,并且在12345678.f
和12345679.f
之间没有单精度浮点数这两个。
(有点)好消息
然而,精确的插值可以表示为 5 个浮点值(在二维情况下为向量)的总和:一个用于公式的结果,一个用于每个操作中的错误 [1] 和一个用于将误差乘以 t
。我不确定这是否对您有用。下面是单精度算法的一维 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)
returns 12345678.f
和 err1
、err2
、err3
和 [= 的结果29=]分别是0.0f
、0.0f
、0.300000011920928955078125f
和0.0f
。事实上,正确的结果是 12345678.300000011920928955078125 不能表示为单精度浮点数。
一个更复杂的例子:exact_lerp(0.23456789553165435791015625f, 7.345678806304931640625f, 0.300000011920928955078125f, &err1, &err2, &err3, &err4)
returns 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] 格雷拉特,斯特夫 (2007)。 Accurate Floating Point Product and Exponentiation.
- [2] Enabling strict floating point mode in GCC
- [3]Semantics of Floating Point Math in GCC