如果在整个过程中使用IEEE-754二进制浮点数的一种格式(例如,基本的32位二进制浮点数,这是C++中常用的
float
格式),并且所有C++运算符以直接简单的方式映射到IEEE-754操作,则
lerp_symmetric(alpha, x0, x1)
(以下简称为
A
)等于
lerp_symmetric(1-alpha, x1, x0)
(
B
)
证明:
如果假设的范围[0, 1]中的alpha大于等于½,则根据Sterbenz引理,1-alpha是精确的(即计算得到的浮点结果等于数学结果,没有舍入误差)。因此,在计算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的w1)是精确的,因为其数学结果是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. 我认为这是一个合理的方法。没有进一步的思考,我不会断言没有可以改进的地方。
2. 不,你不能信任你的编译器。
C++允许在计算浮点数运算时使用多余的精度。因此,即使所有操作数都是float类型,
w0*x0 + w1*x1
可能会使用double、long double或其他精度进行计算。
C++允许缩写,除非禁用,所以
w0*x0 + w1*x1
可以被计算为
fmaf(w0, x0, w1*x1)
,从而对其中一个乘法使用精确算术,但对另一个乘法不使用精确算术。
您可以通过使用以下方法部分解决这个问题:
float w0 = 1.0f - alpha;
float w1 = 1.0f - w0;
float t0 = w0*x0;
float t1 = w1*x1;
return t0+t1;
C++标准要求在赋值和类型转换中丢弃多余的精度。这也适用于函数返回值。(我根据记忆报告了这个和其他C++规范;应该检查标准文件。)因此,即使最初使用了额外的精度,上述每个操作都会将其结果四舍五入为“float”。这将防止缩减。
(通过包含并插入预处理器指令#pragma STDC FP_CONTRACT OFF,还可以禁用缩减。某些编译器可能不支持该功能。)
上述解决方法的一个问题是,值首先被舍入到评估精度,然后再舍入为float
。对于某些数学值,例如某个值x,先将x舍入到double
(或其他精度),然后再舍入为float
与直接将x舍入到float
产生的结果不同。Samuel A. Figueroa del Cid的论文“一个严格支持高级编程语言中IEEE浮点算术标准的框架”指出,在IEEE-754基本64位浮点数(通常用于double
)执行单个乘法或加法操作,然后舍入为32位格式永远不会出现双舍入误差(因为这些操作针对32位格式的输入元素,永远不会生成上述麻烦的x值之一)。1
如果我从记忆中对C++规范的报告正确,那么上述描述的解决方法应该是完整的,只要C++实现在评估浮点表达式时使用名义格式或足够宽的格式来满足Figueroa del Cid给出的要求。
脚注
1 根据Figueroa del Cid,如果x
和y
具有p位有效数字,并且精确计算x+y
或x*y
,然后四舍五入到q位小数,如果将结果再四舍五入到p位小数,将得到与直接将结果四舍五入到p位小数相同的答案,当且仅当p ≤ (q − 1)/2时。这适用于IEEE-754基本32位二进制浮点数(p = 24)和64位(q = 53)。这些格式通常用于float
和double
,在使用它们的C++实现中,上述解决方法应该足够。如果C++实现使用的精度不满足Figueroa del Cid给出的条件,则可能会发生双重舍入错误。
float
精度的IEEE-754算术,我相信您的附加属性就可以得到满足。C和C++允许灵活地评估w0*x0 + w1*x1
。您可以通过编写t0 = w0*x0; t1 = w1*x1; return t0+t1;
来有些抵消这种情况,因为赋值必须丢弃多余的精度。然后,语言规则仍然允许计算例如t0 = x0*x0
的多余精度,然后缩小为float
。但是,我相信有一篇论文表明,只要多余精度足够大,这等效于在float
中进行计算。 - Eric Postpischil