c - 从 double 参数开始的80位扩展精度计算的属性

这是插值函数的两种实现。参数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);
}

在具有80位long double的严格IEEE 754平台上，interpol_64()中的所有计算均根据IEEE 754 double 进行，而interpol_80()中的所有计算均以80位扩展精度进行。
程序打印:

u2: 0x1.8f5c28f5c28f6p-3
y64:0x1.8f5c28f5c28f5p-3
y80:0x1.8f5c28f5c28f6p-3

我对属性“感兴趣”(该函数返回的结果始终在u2和u3之间)。如上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的80位。

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
最后要研究的情况是u2 因此，我们所做的断言不再成立:

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

优质教育

为了完整起见，我们应该注意不正常的输入(逐渐下溢)，但是我希望您在压力测试中不会那么恶毒。我不会证明这些发生了什么...

编辑:

这是后续操作，因为以下断言有点近似，并且在0

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

我没有证明原始要求，但还没有到此为止。

关于c - 从 double 参数开始的80位扩展精度计算的属性，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/13725802/

U2

c - 从 double 参数开始的80位扩展精度计算的属性