下面是插值功能的两种实现方法。 参数u1之间总是0.和1.
#include <stdio.h>
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);
}
在严格的IEEE 754与平台80位long double S，在所有计算interpol_64()是根据IEEE 754双精度完成，并且在interpol_80()在80位的扩展精度。 该程序打印：
u2: 0x1.8f5c28f5c28f6p-3
y64:0x1.8f5c28f5c28f5p-3
y80:0x1.8f5c28f5c28f6p-3
我感兴趣的财产“由该函数返回的结果总是在两者之间u2和u3 ”。 此属性为false的interpol_64()如在所示的值main()的上方。
是否拥有房产的机会成为真正的interpol_80() 如果不是，什么是一个反例？ 它是否有助于如果我们知道u2 != u3或它们之间的最小距离？ 是否有一个来确定用于在该属性将被保证是真中间计算的有效数宽度的方法？
编辑：对所有我试过的随机值，当中间计算在扩展精度都在内部完成的物业持有。 如果interpol_80()花了long double参数，这将是比较容易建立一个反例，但这里的问题是，特别是约接受一个函数double参数。 这使得它更难建立一个反例，如果有一个。
注：编译器产生的x87指令可能产生相同的代码interpol_64()和interpol_80()但是这是切我的问题。
是的，interpol_80（）是安全的，让我们来示范一下。
这个问题指出，输入是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
转换为U1，U2，U3的80bits是准确的了。
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是关闭，因此刑警-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
是完整的，我们应该关心的非正规输入（逐渐下溢），但我希望你不会有压力测试是恶性。 我不会表现出那些发生了什么？
编辑 ：
这里是一个后续如以下断言有点近似和产生一些评论当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
我没有证明原件要求，但不会太远......
在失精度的主要来源interpol_64是乘法。 乘以2 53位的尾数产生一个105-或106位（取决于高比特是否携带）尾数。 这是太大，无法在80位扩展精度值，所以在一般情况下，你还必须在80位版本的失精度。 准确量化当它发生是非常困难的; 这是容易说的最多的是，当舍入误差积累它发生。 请注意，增加这两个词的时候也是一个小舍入一步。
大多数人都可能只是解决这个问题，就像一个功能：
double interpol_64(double u1, double u2, double u3)
{
return u2 + u1 * (u3 - u2);
}
但它看起来像你要找的洞察四舍五入问题，而不是一个更好的实现。
文章来源: Properties of 80-bit extended precision computations starting from double precision arguments
