这不能回答最初提出的问题，但如果您是处理类似问题的程序员，这个答案可能会对您有所帮助。
我真的不明白感知到的困难在哪里。提供严格的 IEEE-754 binary64 语义，同时仅限于 80387 浮点数学，并保留 80 位长双精度计算，似乎遵循 GCC-4.6.3 和 clang-3.0 的明确规定的 C99 转换规则（基于 LLVM 3.0）。
编辑添加：然而，Pascal Cuoq 是正确的：gcc-4.6.3 或 clang-llvm-3.0 实际上都没有为 '387 浮点数学正确执行这些规则。给定正确的编译器选项，规则会正确应用于编译时评估的表达式，但不适用于运行时表达式。有一些解决方法，在下面的休息后列出。
我编写分子动力学模拟代码，并且非常熟悉可重复性/可预测性要求，并且希望尽可能保持最大精度，所以我确实声称我知道我在说什么。这个答案应该表明这些工具存在并且易于使用；问题源于不了解或不使用这些工具。
（我喜欢的一个首选示例是 Kahan 求和算法。使用 C99 和适当的强制转换（将强制转换添加到例如 Wikipedia 示例代码），根本不需要任何技巧或额外的临时变量。无论编译器优化级别如何，实现都可以工作，包括-O3和-Ofast。）
C99 明确声明（在例如 5.4.2.2 中）强制转换和赋值都删除了所有额外的范围和精度。这意味着您可以long double通过将计算期间使用的临时变量定义为 来使用算术long double，也可以将输入变量转换为该类型；每当需要 IEEE-754 binary64 时，只需转换为double.
在 '387 上，强制转换在上述两个编译器上生成一个赋值和一个负载；这确实将 80 位值正确舍入为 IEEE-754 binary64。在我看来，这个成本是非常合理的。所花费的确切时间取决于架构和周围的代码；通常它可以并且可以与其他代码交错以将成本降低到可以忽略的水平。当 MMX、SSE 或 AVX 可用时，它们的寄存器与 80 位 80387 寄存器是分开的，并且通常通过将值移动到 MMX/SSE/AVX 寄存器来完成转换。
（我更喜欢生产代码对临时变量使用特定的浮点类型，比如说tempdouble或这样，以便可以将其定义为double或long double取决于所需的架构和速度/精度权衡。）
简而言之：
不要仅仅因为所有变量和文字常量都是精确的(expression)。double把它写成(double)(expression)你想要double精确的结果。
这也适用于复合表达式，并且有时可能会导致具有多个强制转换级别的笨拙表达式。
如果您有expr1并且expr2希望以 80 位精度计算，但还需要先将每个四舍五入到 64 位的乘积，请使用
long double expr1;
long double expr2;
double product = (double)(expr1) * (double)(expr2);
注意，product计算为两个 64 位值的乘积；未以 80 位精度计算，然后向下舍入。以 80 位精度计算乘积，然后向下舍入，将是
double other = expr1 * expr2;
或者，添加描述性的演员来告诉你到底发生了什么，
double other = (double)((long double)(expr1) * (long double)(expr2));
这应该是显而易见的，product而且other经常有所不同。
如果您确实使用混合的 32 位/64 位/80 位/128 位浮点值，C99 转换规则只是您必须学会使用的另一个工具。float确实，如果您混合使用 binary32 和 binary64 浮点数（以及double在大多数架构上），您会遇到完全相同的问题！
也许重写 Pascal Cuoq 的探索代码，以正确应用强制转换规则，使这更清楚？
#include <stdio.h>
#define TEST(eq) printf("%-56s%s\n", "" # eq ":", (eq) ? "true" : "false")
int main(void)
{
double d = 1.0 / 10.0;
long double ld = 1.0L / 10.0L;
printf("sizeof (double) = %d\n", (int)sizeof (double));
printf("sizeof (long double) == %d\n", (int)sizeof (long double));
printf("\nExpect true:\n");
TEST(d == (double)(0.1));
TEST(ld == (long double)(0.1L));
TEST(d == (double)(1.0 / 10.0));
TEST(ld == (long double)(1.0L / 10.0L));
TEST(d == (double)(ld));
TEST((double)(1.0L/10.0L) == (double)(0.1));
TEST((long double)(1.0L/10.0L) == (long double)(0.1L));
printf("\nExpect false:\n");
TEST(d == ld);
TEST((long double)(d) == ld);
TEST(d == 0.1L);
TEST(ld == 0.1);
TEST(d == (long double)(1.0L / 10.0L));
TEST(ld == (double)(1.0L / 10.0));
return 0;
}
GCC 和 clang 的输出是
sizeof (double) = 8
sizeof (long double) == 12
Expect true:
d == (double)(0.1): true
ld == (long double)(0.1L): true
d == (double)(1.0 / 10.0): true
ld == (long double)(1.0L / 10.0L): true
d == (double)(ld): true
(double)(1.0L/10.0L) == (double)(0.1): true
(long double)(1.0L/10.0L) == (long double)(0.1L): true
Expect false:
d == ld: false
(long double)(d) == ld: false
d == 0.1L: false
ld == 0.1: false
d == (long double)(1.0L / 10.0L): false
ld == (double)(1.0L / 10.0): false
除了最近版本的 GCC 将右手边提升ld == 0.1为 long double first（即 to ld == 0.1L），yieldingtrue和 SSE/AVX的右手边long double是 128 位的。
对于纯 '387 测试，我使用了
gcc -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
clang -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
具有各种优化标志组合...，包括-fomit-frame-pointer、-O0、-O1、-O2、-O3和-Os。
使用任何其他标志或 C99 编译器应该会导致相同的结果，除了long double大小（以及ld == 1.0当前的 GCC 版本）。如果您遇到任何差异，我将非常感谢您了解它们；我可能需要警告我的用户使用此类编译器（编译器版本）。请注意，Microsoft 不支持 C99，因此我对它们完全不感兴趣。
Pascal Cuoq 确实在下面的评论链中提出了一个有趣的问题，我没有立即意识到。
在计算表达式时，GCC 和 clang with-mfpmath=387都指定使用 80 位精度计算所有表达式。这导致例如
7491907632491941888 = 0x1.9fe2693112e14p+62 = 110011111111000100110100100110001000100101110000101000000000000
5698883734965350400 = 0x1.3c5a02407b71cp+62 = 100111100010110100000001001000000011110110111000111000000000000
7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000
产生不正确的结果，因为二进制结果中间的一串 1 正好在 53 位和 64 位尾数之间（分别为 64 位和 80 位浮点数）。所以，虽然预期的结果是
42695510550671088819251326462451515392 = 0x1.00f6d98d0a42fp+125 = 100000000111101101101100110001101000010100100001011110000000000000000000000000000000000000000000000000000000000000000000000000
用just得到的结果-std=c99 -m32 -mno-sse -mfpmath=387是
42695510550671098263984292201741942784 = 0x1.00f6d98d0a43p+125 = 100000000111101101101100110001101000010100100001100000000000000000000000000000000000000000000000000000000000000000000000000000
理论上，您应该能够通过使用选项告诉 gcc 和 clang 强制执行正确的 C99 舍入规则
-std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard
但是，这只影响编译器优化的表达式，并且似乎根本无法修复 387 处理。如果您将 egclang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test与Pascal Cuoq 的示例程序test.c一起使用，您将根据 IEEE-754 规则获得正确的结果——但这只是因为编译器优化了表达式，根本不使用 387。
简单地说，而不是计算
(double)d1 * (double)d2
gcc 和 clang 实际上都告诉 '387 进行计算
(double)((long double)d1 * (long double)d2)
这确实是我相信这是一个影响 gcc-4.6.3 和 clang-llvm-3.0 的编译器错误，并且很容易复制。（Pascal Cuoq 指出这FLT_EVAL_METHOD=2意味着对双精度参数的操作总是以扩展精度完成，但我看不出任何合理的理由——除了必须重写libm'387 的部分内容——在 C99 中这样做并考虑到 IEEE- 754 规则可以通过硬件实现！毕竟，正确的操作很容易被编译器实现，通过修改 '387 控制字以匹配表达式的精度。并且，考虑到应该强制这种行为的编译器选项-std=c99 -ffloat-store -fexcess-precision=standard—— - 如果没有意义FLT_EVAL_METHOD=2行为实际上是需要的，也不存在向后兼容性问题。）重要的是要注意，给定正确的编译器标志，在编译时评估的表达式确实会得到正确的评估，并且只有在运行时评估的表达式才会得到不正确的结果。
最简单的解决方法和可移植的解决方法是使用fesetround(FE_TOWARDZERO)(from fenv.h) 将所有结果四舍五入为零。
在某些情况下，向零舍入可能有助于可预测性和病理情况。特别是，对于像 的区间x = [0,1)，向零舍入意味着通过舍入永远不会达到上限；如果您评估例如分段样条曲线，这很重要。
对于其他舍入模式，您需要直接控制 387 硬件。
您可以使用__FPU_SETCW()from#include <fpu_control.h>或打开代码。例如precision.c：
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#define FP387_NEAREST 0x0000
#define FP387_ZERO 0x0C00
#define FP387_UP 0x0800
#define FP387_DOWN 0x0400
#define FP387_SINGLE 0x0000
#define FP387_DOUBLE 0x0200
#define FP387_EXTENDED 0x0300
static inline void fp387(const unsigned short control)
{
unsigned short cw = (control & 0x0F00) | 0x007f;
__asm__ volatile ("fldcw %0" : : "m" (*&cw));
}
const char *bits(const double value)
{
const unsigned char *const data = (const unsigned char *)&value;
static char buffer[CHAR_BIT * sizeof value + 1];
char *p = buffer;
size_t i = CHAR_BIT * sizeof value;
while (i-->0)
*(p++) = '0' + !!(data[i / CHAR_BIT] & (1U << (i % CHAR_BIT)));
*p = '\0';
return (const char *)buffer;
}
int main(int argc, char *argv[])
{
double d1, d2;
char dummy;
if (argc != 3) {
fprintf(stderr, "\nUsage: %s 7491907632491941888 5698883734965350400\n\n", argv[0]);
return EXIT_FAILURE;
}
if (sscanf(argv[1], " %lf %c", &d1, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[1]);
return EXIT_FAILURE;
}
if (sscanf(argv[2], " %lf %c", &d2, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[2]);
return EXIT_FAILURE;
}
printf("%s:\td1 = %.0f\n\t %s in binary\n", argv[1], d1, bits(d1));
printf("%s:\td2 = %.0f\n\t %s in binary\n", argv[2], d2, bits(d2));
printf("\nDefaults:\n");
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to nearest integer:\n");
fp387(FP387_EXTENDED | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to nearest integer:\n");
fp387(FP387_DOUBLE | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to zero:\n");
fp387(FP387_EXTENDED | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to zero:\n");
fp387(FP387_DOUBLE | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
return 0;
}
使用clang-llvm-3.0编译运行，得到正确结果，
clang -std=c99 -m32 -mno-sse -mfpmath=387 -O3 -W -Wall precision.c -o precision
./precision 7491907632491941888 5698883734965350400
7491907632491941888: d1 = 7491907632491941888
0100001111011001111111100010011010010011000100010010111000010100 in binary
5698883734965350400: d2 = 5698883734965350400
0100001111010011110001011010000000100100000001111011011100011100 in binary
Defaults:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Extended precision, rounding to nearest integer:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Double precision, rounding to nearest integer:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Extended precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Double precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
换句话说，您可以通过fp387()设置精度和舍入模式来解决编译器问题。
缺点是某些数学库（libm.a, libm.so）可能会在假设中间结果始终以 80 位精度计算的情况下编写。至少fpu_control.hx86_64 上的 GNU C 库有注释"libm requires extended precision"。幸运的是，您可以从例如 GNU C 库中获取 '387 实现，并在头文件中实现它们或编写 known-to-work libm（如果您需要该math.h功能）；事实上，我想我可能会在那里提供帮助。