【发布时间】：2012-11-23 10:19:58
【问题描述】：
这里有两种插值函数的实现。参数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 doubles 的严格 IEEE 754 平台上，interpol_64() 中的所有计算均根据 IEEE 754 双精度进行，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() 生成相同的代码，但这与我的问题无关。
【问题讨论】：
-
你确定这个程序真的使用 80 位精度吗？ IIRC 现代 Intel / AMD 机器内置 128 fp 单元，随 SSE 和朋友一起提供。
-
@FUZxxl “128 位 FP 单元”表示两个双精度或 4 个单精度数字的向量。但要回答你的问题，是的，我敢肯定。大会在这里：pastebin.com/GaM20WZS
-
内容和演示都+1
标签： c floating-point ieee-754 extended-precision