是否有文档描述Clang如何处理过多的浮点精度？

当允许使用的唯一的浮点指令是387时，几乎不可能（*）以合理的成本提供严格的IEEE 754语义。当希望保持FPU工作在完整的64位有效位上时，这是特别困难的，所以long double精度型可用于扩展精度。通常的“解决scheme”是以唯一可用的精度进行中间计算，并在或多或less明确定义的情况下将其转换为较低的精度。

最近的GCC版本根据Joseph S. Myers在2008年发布的GCC邮件列表中的解释，处理了中间计算中的超额精度。就我所知，这个描述使得一个程序编译成gcc -std=c99 -mno-sse2 -mfpmath=387完全可以预测的。如果偶然的话，它不是，这是一个错误，它将是固定的：Joseph S. Myers在他的文章中表示的意图是使其可预测。

是否logging了Clang如何处理过度的精度（比如使用选项-mno-sse2 ）以及在哪里？

（*）编辑：这是一个夸张。当允许将x87 FPUconfiguration为使用53位有效位时，它会稍微令人烦恼，但是仿真二进制64 并不困难。

在下面的R ..的评论之后，下面是我和Clang最近版本的简短交互的日志：

 Hexa:~ $ clang -v Apple clang version 4.1 (tags/Apple/clang-421.11.66) (based on LLVM 3.1svn) Target: x86_64-apple-darwin12.4.0 Thread model: posix Hexa:~ $ cat fem.c #include <stdio.h> #include <math.h> #include <float.h> #include <fenv.h> double x; double y = 2.0; double z = 1.0; int main(){ x = y + z; printf("%d\n", (int) FLT_EVAL_METHOD); } Hexa:~ $ clang -std=c99 -mno-sse2 fem.c Hexa:~ $ ./a.out 0 Hexa:~ $ clang -std=c99 -mno-sse2 -S fem.c Hexa:~ $ cat fem.s … movl $0, %esi fldl _y(%rip) fldl _z(%rip) faddp %st(1) movq _x@GOTPCREL(%rip), %rax fstpl (%rax) …

这不能回答最初提出的问题，但如果你是一个处理类似问题的程序员，这个答案可能会帮助你。

我真的没有看到困难在哪里。提供严格的IEEE-754二进制64语义，同时限制为80387浮点math，保留80位长的双重计算，似乎在GCC-4.6.3和clang-3.0（基于LLVM 3.0）。

编辑补充：然而，Pascal Cuoq是正确的：gcc-4.6.3或clang-llvm-3.0实际上对387浮点math实际执行这些规则。给定正确的编译器选项，规则将正确应用于编译时评估的expression式，但不适用于运行时expression式。有一些解决方法，在下面的分支之后列出。

我做分子动力学模拟的代码，对可重复性/可预测性的要求非常熟悉，也希望尽可能保持最大的精度，所以我确实声称我知道我在这里谈论什么。这个答案应该表明这些工具存在并且易于使用; 这些问题是由于不知道或不使用这些工具而产生的。

（我喜欢的一个首选例子是Kahan求和algorithm，通过C99和适当的转换（例如添加到维基百科的示例代码），根本不需要任何技巧或额外的临时variables，不pipe编译器优化级别如何， -O3和-Ofast 。）

C99明确指出（例如在5.4.2.2中）铸造和分配都会消除所有额外的范围和精度。这意味着你可以使用long double算术，在计算过程中使用临时variables为long double ，也可以将inputvariables转换为该types; 每当需要IEEE-754二进制64时，只需转换成双精度。

在'387上，演员们在上面的编译器上产生了一个任务和一个负载。这确实将80位值转换为IEEE-754二进制64。我认为这个成本是非常合理的。确切的时间取决于体系结构和周围的代码; 通常情况下，它可以与其他代码交错，从而将成本降至可以忽略的程度。当MMX，SSE或AVX可用时，它们的寄存器与80位80387寄存器分开，通常通过将值移到MMX / SSE / AVX寄存器来完成。

（对于临时variables，我更喜欢使用特定的浮点types（比如tempdouble或者这样的types）的生产代码，以便可以根据所需的体系结构和速度/精度权衡来定义为double或long double 。）

简而言之：

不要假设(expression)是double精度的，因为所有的variables和文字常量都是。把它写成(double)(expression)如果你想得到double精度的结果。

这也适用于复合expression式，并且有时可能会导致具有不同级别的强制转换的笨拙的expression。

如果你想要以80位精度计算expr1和expr2 ，而且还需要每个舍入到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转换规则只是您必须学习使用的另一个工具。真的，如果你混合使用binary32和binary64浮点数（在大多数架构上使用float和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和叮当的输出是

 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的第一个位置（即ld == 0.1L ），这个结果是true ，而对于SSE / AVX， long double是128位的。

对于纯粹的'387testing，我用了

 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大小（当前GCC版本和ld == 1.0 ）。如果您遇到任何分歧，我会非常感激听到他们的消息。我可能需要警告这些编译器（编译器版本）的用户。请注意，微软不支持C99，所以他们对我完全不感兴趣。

Pascal Cuoq在下面的评论链中提出了一个有趣的问题，我没有马上意识到。

在评估expression式时，GCC和clang与-mfpmath=387指定所有expression式都使用80位精度进行计算。这导致例如

 7491907632491941888 = 0x1.9fe2693112e14p+62 = 110011111111000100110100100110001000100101110000101000000000000 5698883734965350400 = 0x1.3c5a02407b71cp+62 = 100111100010110100000001001000000011110110111000111000000000000 7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000

会产生不正确的结果，因为在二进制结果中间的那串string只是53和64位尾数（分别为64和80位浮点数）之间的差异。所以，虽然预期的结果是

 42695510550671088819251326462451515392 = 0x1.00f6d98d0a42fp+125 = 100000000111101101101100110001101000010100100001011110000000000000000000000000000000000000000000000000000000000000000000000000

用-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

但是，这只会影响编译器优化的expression式，而且似乎不能修复387处理。如果你使用eg clang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test用test.c作为Pascal Cuoq的例子程序，您将得到每个IEEE-754规则的正确结果 – 但仅仅是因为编译器优化了expression式，而不是使用387。

简单地说，而不是计算

 (double)d1 * (double)d2

海湾合作委员会和铿锵实际上告诉'387计算

 (double)((long double)d1 * (long double)d2)

~~这确实是~~ 我相信这是一个影响gcc-4.6.3和clang-llvm-3.0的编译器bug，也是一个很容易被复制的bug。（Pascal Cuoq指出， FLT_EVAL_METHOD=2意味着对双精度参数的操作总是以FLT_EVAL_METHOD=2的精度完成的，但是我看不到任何理智的原因 – 除了必须在387上重写libm部分外 – 在C99中这样做并考虑到IEEE-754规则是可以通过硬件来实现的！毕竟，编译器很容易实现正确的操作，通过修改'387控制字来匹配expression式的精度，而且，行为 – -std=c99 -ffloat-store -fexcess-precision=standard – 如果FLT_EVAL_METHOD=2行为实际上是需要的，那么也没有意义，也不存在向后兼容性问题。）重要的是要注意，编译器标志，编译时评估的expression式得到正确的评估，并且只有在运行时评估的expression式得到不正确的结果。

最简单的解决方法是使用fesetround(FE_TOWARDZERO) （来自fenv.h ）将所有结果四舍五入到零。

在某些情况下，趋于零可能有助于预测性和病理情况。特别是，对于x = [0,1)这样的区间，舍入为零意味着通过舍入决不会达到上限; 重要的是如果你评估例如分段样条。

对于其他舍入模式，您需要直接控制387硬件。

您可以使用#include <fpu_control.h> __FPU_SETCW() ，或者打开它的代码。例如， 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()来设置精度和舍入模式来解决编译器问题。

不利的一面是，一些math库（ libm.a ， libm.so ）可能会被假设为中间结果总是以80位精度进行计算。至lessx86_64上的GNU C库fpu_control.h有注释“libm需要扩展精度” 。幸运的是，如果您需要math.hfunction，您可以从例如GNU C库中获取'387实现，并将其实现为头文件或写入已知工作的libm 。实际上，我想我可以在那里帮忙。

为了logging，下面是我通过实验发现的。以下程序显示了使用Clang编译时的各种行为：

 #include <stdio.h> int r1, r2, r3, r4, r5, r6, r7; double ten = 10.0; int main(int c, char **v) { r1 = 0.1 == (1.0 / ten); r2 = 0.1 == (1.0 / 10.0); r3 = 0.1 == (double) (1.0 / ten); r4 = 0.1 == (double) (1.0 / 10.0); ten = 10.0; r5 = 0.1 == (1.0 / ten); r6 = 0.1 == (double) (1.0 / ten); r7 = ((double) 0.1) == (1.0 / 10.0); printf("r1=%d r2=%d r3=%d r4=%d r5=%d r6=%d r7=%d\n", r1, r2, r3, r4, r5, r6, r7); }

结果因优化级别而异：

 $ clang -v Apple LLVM version 4.2 (clang-425.0.24) (based on LLVM 3.2svn) $ clang -mno-sse2 -std=c99 tc && ./a.out r1=0 r2=1 r3=0 r4=1 r5=1 r6=0 r7=1 $ clang -mno-sse2 -std=c99 -O2 tc && ./a.out r1=0 r2=1 r3=0 r4=1 r5=1 r6=1 r7=1

在-O2上区分r5和r6的cast (double)对-O0和variablesr3和r4没有影响。 r1与r5在所有优化级别上都是不同的，而r6只与r6不同。