从双精度参数开始的80位扩展精度计算的属性
以下是插值函数的两种实现。 参数u1
始终在0.
和1.
之间。
#include double interpol_64(double u1, double u2, double u3) { return u2 * (1.0 - u1) + u1 * u3; } double interpol_80(double u1, double u2, double u3) { return u2 * (1.0 - (long double)u1) + u1 * (long double)u3; } int main() { double y64,y80,u1,u2,u3; u1 = 0.025; u2 = 0.195; u3 = 0.195; y64 = interpol_64(u1, u2, u3); y80 = interpol_80(u1, u2, u3); printf("u2: %a\ny64:%a\ny80:%a\n", u2, y64, y80); }
在具有80位long double
s的严格IEEE 754平台上, interpol_64()
中的所有计算均根据IEEE 754双精度进行,而interpol_64()
则以80位扩展精度进行。 程序打印:
u2: 0x1.8f5c28f5c28f6p-3 y64:0x1.8f5c28f5c28f5p-3 y80:0x1.8f5c28f5c28f6p-3
我对属性感兴趣“函数返回的结果总是介于u2
和u3
之间”。 此属性为interpol_64()
false,如上面main()
的值所示。
该物业是否有机会成为interpol_80()
真实情况? 如果不是,反例是什么? 如果我们知道u2 != u3
或者它们之间有最小距离会有帮助吗? 是否有一种方法可以确定中间计算的有效位宽度,在该计算中,属性将保证为真?
编辑:在我尝试的所有随机值上,当在内部以扩展精度完成中间计算时,属性保持不变。 如果interpol_80()
采用long double
参数,那么构建一个反例也相对容易,但这里的问题具体是关于一个带double
参数的函数。 这使得构建反例(如果有的话)变得更加困难。
注意:生成x87指令的编译器可能会为interpol_64()
和interpol_80()
生成相同的代码,但这与我的问题相关。
是的,interpol_80()是安全的,让我们来certificate一下。
问题表明输入是64位浮点数
rnd64(ui) = ui
结果是完全的(假设*和+是数学运算)
r = u2*(1-u1)+(u1*u3)
舍入为64位浮点数的最佳返回值为
r64 = rnd64(r)
因为我们有这些属性
u2 <= r <= u3
这是有保证的
rnd64(u2) <= rnd64(r) <= rnd64(u3) u2 <= r64 <= u3
转换为80位的u1,u2,u3也是准确的。
rnd80(ui)=ui
现在,假设0 <= u2 <= u3
,然后执行不精确的浮点运算导致最多4个舍入错误:
rf = rnd(rnd(u2*rnd(1-u1)) + rnd(u1*u3))
假设舍入到最接近的偶数,则精确值最多为2 ULP。 如果使用64位浮点数或80位浮点数执行舍入:
r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r) r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r)
rf64
可以关闭2 ulp所以interpol-64()是不安全的,但rnd64( rf80 )
怎么样?
我们可以说:
rnd64(r - 2 ulp80(r)) <= rnd64(rf80) <= rnd64(r + 2 ulp80(r))
因为0 <= u2 <= u3
,那么
ulp80(u2) <= ulp80(r) <= ulp80(r3) rnd64(u2 - 2 ulp80(u2)) <= rnd64(r - 2 ulp80(r)) <= rnd64(rf80) rnd64(u3 + 2 ulp80(u3)) >= rnd64(r + 2 ulp80(r)) >= rnd64(rf80)
幸运的是,就像范围内的每个数字(u2-ulp64(u2)/2 , u2+ulp64(u2)/2)
我们得到的
rnd64(u2 - 2 ulp80(u2)) = u2 rnd64(u3 + 2 ulp80(u3)) = u3
因为ulp80(x)=ulp62(x)/2^(64-53)
因此,我们得到证据
u2 <= rnd64(rf80) <= u3
对于u2 <= u3 <= 0,我们可以轻松应用相同的证明。
最后要研究的案例是u2 <= 0 <= u3。 如果我们减去2个大值,那么结果可以达到ulp(大)/ 2关闭而不是ulp(大 - 大)/ 2 ......
因此,我们所做的断言不再适用:
r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r)
幸运的是, u2 <= u2*(1-u1) <= 0 <= u1*u3 <= u3
并且在舍入后保留
u2 <= rnd(u2*rnd(1-u1)) <= 0 <= rnd(u1*u3) <= u3
因此,由于增加的数量符号相反:
u2 <= rnd(u2*rnd(1-u1)) + rnd(u1*u3) <= u3
四舍五入之后也是如此,所以我们可以再次保证
u2 <= rnd64( rf80 ) <= u3
QED
为了完整,我们应该关注非正规输入(逐渐下溢),但我希望你不会对压力测试产生恶意。 我不会certificate那些会发生什么......
编辑 :
这是一个后续行动,因为以下断言有点近似,并在0 <= u2 <= u3时生成一些注释
r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r)
我们可以写出以下不等式:
rnd(1-u1) <= 1 rnd(1-u1) <= 1-u1+ulp(1)/4 u2*rnd(1-u1) <= u2 <= r u2*rnd(1-u1) <= u2*(1-u1)+u2*ulp(1)/4 u2*ulp(1) < 2*ulp(u2) <= 2*ulp(r) u2*rnd(1-u1) < u2*(1-u1)+ulp(r)/2
对于下一轮操作,我们使用
ulp(u2*rnd(1-u1)) <= ulp(r) rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(u2*rnd(1-u1))/2 rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(r)/2 rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)
对于总和的第二部分,我们有:
u1*u3 <= r rnd(u1*u3) <= u1*u3 + ulp(u1*u3)/2 rnd(u1*u3) <= u1*u3 + ulp(r)/2 rnd(u2*rnd(1-u1))+rnd(u1*u3) < u2*(1-u1)+u1*u3 + 3*ulp(r)/2 rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 3*ulp(r)/2 + ulp(r+3*ulp(r)/2)/2 ulp(r+3*ulp(r)/2) <= 2*ulp(r) rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 5*ulp(r)/2
我没有certificate原来的说法,但不是那么远......
interpol_64
精度损失的主要来源是乘法。 乘以两个53位尾数产生105或106位(取决于高位是否携带)尾数。 这太大而不适合80位扩展精度值,因此通常,80位版本中也会出现精度损失。 准确地量化它何时发生是非常困难的; 最容易说的是,当舍入错误累积时会发生这种情况。 请注意,添加这两个术语时,还会有一个小的舍入步骤。
大多数人可能会用以下函数解决这个问题:
double interpol_64(double u1, double u2, double u3) { return u2 + u1 * (u3 - u2); }
但看起来您正在寻找对四舍五入问题的洞察力,而不是更好的实施。