我有一个功能:
float lerp(float alpha, float x0, float x1) {
return (1.0f - alpha) * x0 + alpha * x1;
}
对于那些没有见过它的人来说,这比x0 + (x1-x0) * alpha
更好,因为后者并不保证lerp(1.0f, x0, x1) == x1
。
现在,我希望我的lerp
函数有一个额外的属性:我喜欢lerp(alpha, x0, x1) == lerp(1-alpha, x1, x0)
。(至于为什么:这是一个更复杂的函数的玩具示例。)我想出的似乎有效的解决方案是
float lerp_symmetric(float alpha, float x0, float x1) {
float w0 = 1.0f - alpha;
float w1 = 1.0f - w0;
return w0 * x0 + w1 * x1;
}
这种双减法的效果是在零和一附近舍入,所以如果alpha = std::nextafter(0)
(1.4012985e-45),那么1 - alpha == 1
,所以1 - (1-alpha) == 0
。据我所知,1.0f - x == 1.0f - (1.0f - (1.0f - x))
总是正确的。它似乎也有w0 + w1 == 1.0f
的效果。
问题:
1.这做法是否合理呢?
1.我可以信任我的编译器做我想做的事吗?特别是,我知道在Windows上它有时会对部分结果使用更高的精度,我知道编译器被允许做一些代数运算;显然1-(1-x)==x代数。
这是在C++11中使用Clang,VisualStudio和gcc。
1条答案
按热度按时间btqmn9zl1#
如果始终使用IEEE-754二进制浮点的一种格式(例如,基本的32位二进制,通常用于C++
float
的格式),并且所有C++运算符都以直接和简单的方式Map到IEEE-754操作,则lerp_symmetric(alpha, x0, x1)
(以下称为A
)等于lerp_symmetric(1-alpha, x1, x0)
(B
)证明:
alpha
,我们假设它在[0,1]中,大于或等于½,则1-alpha
通过Sterbenz'引理是精确的。(通过“精确”,我们的意思是计算的浮点结果等于数学结果;因此,在计算A
时,w0
是精确的,因为它是1-alpha
,w1
是精确的,因为它的数学结果是alpha
,所以它是精确可表示的。在计算B
时,w0
是精确的,因为它的数学结果是alpha
,w1
是精确的,因为它又是1-alpha
。alpha
小于½,则1-alpha
可能有一些舍入误差。假设结果为beta
。然后,在A
中,w0
是beta
。现在½ ≤beta
,所以Sterbenz引理适用于w1 = 1.0f - w0
的求值,所以w1
是精确的(并且等于1-beta
的数学结果)。而且,在B
中,w0
是精确的,再次通过Sterbenz引理,并且等于A
的w1
,并且w1
(B
)是精确的,因为它的数学结果是beta
,这是可以精确表示的。现在我们可以看到
A
中的w0
等于B
中的w1
,A
中的w1
等于B
中的w0
。在上述任何一种情况下,让beta
为1-alpha
,则A
和B
分别返回(1-beta)*x0 + beta*x1
和beta*x1 + (1-beta)*x0
。IEEE-754加法是可交换的(NaN有效载荷除外),因此A
和B
返回相同的结果。回答问题:
1.我认为这是一个合理的方法。如果没有进一步的思考,我不会Assert没有可以改进的地方。
1.不,你不能信任你的编译器:
float
,w0*x0 + w1*x1
也可以使用double
、long double
或其他精度进行计算。w0*x0 + w1*x1
可以被计算为fmaf(w0, x0, w1*x1)
,因此对其中一个乘法使用精确的算术,而不是另一个。您可以使用以下命令部分解决此问题:
C标准要求在赋值和强制转换中丢弃过多的精度。这扩展到函数返回。(我从记忆中报告这个和其他C规范;因此,即使最初使用了额外的精度,上面的每一个都将其结果四舍五入到
float
。这将防止收缩。(One还应该能够通过包含
<cmath>
并插入预处理器指令#pragma STDC FP_CONTRACT OFF
来禁用收缩。有些编译器可能不支持。)上述解决方法的一个问题是,值首先四舍五入到评估精度,然后四舍五入到
float
。对于这样的值 x,首先将 x 四舍五入到double
(或其他精度),然后再四舍五入到float
,会产生与直接将 x 四舍五入到float
不同的结果。本文是SamuelA。Figueroa del Cid证明,在IEEE-754基本64位浮点(通常用于double
)中计算乘法或加法的单个运算,然后舍入到32位格式,永远不会产生双重舍入错误(因为这些运算,给定的输入是32位格式的元素,永远不会产生上述麻烦的 x 值之一)。如果我对从内存中报告的C规范是正确的,那么只要C实现使用标称格式或足够宽的格式来计算浮点表达式以满足Figueroa del Cid给出的要求,上述解决方案就应该是完整的。
脚注
1根据Figueroa del Cid,如果
x
和y
具有 p 位有效数,并且x+y
或x*y
被精确计算,然后四舍五入到 q 位,则第二次四舍五入到 p 位将具有与如果 p ≤(q − 1)/2则结果直接四舍五入到 p 位相同的答案。这对于IEEE-754基本32位二进制浮点(p = 24)和64位(q = 53)是满足的。这些格式通常用于float
和double
,上面的解决方法在使用它们的C实现中应该足够了。如果C实现使用不满足Figueroa del Cid给出的条件的精度来计算float
,则可能会发生双重舍入错误。