11

这是插值函数的两种实现。参数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 双精度和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()and生成相同的代码interpol_80()，但这与我的问题无关。

4

2 回答 2

3

是的，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

转换为 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(big)/2 off 而不是 ulp(big-big)/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

量子点

为了完整起见，我们应该关注非正规输入（逐渐下溢），但我希望你不会对压力测试那么恶毒。我不会演示这些会发生什么...

编辑：

这是一个后续，因为以下断言有点近似，并且在 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

我没有证明最初的说法，但不是那么远......

于 2012-12-05T21:24:55.813 回答
2

精度损失的主要来源interpol_64是乘法。将两个 53 位尾数相乘产生一个 105 位或 106 位（取决于高位是否携带）尾数。这对于 80 位扩展精度值来说太大了，因此通常在 80 位版本中也会出现精度损失。准确量化它何时发生是非常困难的；最容易说的是，当舍入误差累积时会发生这种情况。请注意，添加这两个术语时还有一个小的舍入步骤。

大多数人可能会使用以下功能解决此问题：

double interpol_64(double u1, double u2, double u3)
{ 
  return u2 + u1 * (u3 - u2);
}

但看起来您正在寻找对舍入问题的洞察力，而不是更好的实现。

于 2012-12-05T15:24:05.727 回答