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

Question

当允许使用的唯一浮点指令是387个时,几乎不可能(*)以合理的成本提供严格的IEEE 754语义.当希望保持FPU在完整的64位有效数字上工作时,这一点特别困难,因此该long double类型可用于扩展精度.通常的"解决方案"是以唯一可用的精度进行中间计算,并在或多或少明确定义的场合转换为较低的精度.

根据Joseph S. Myers在2008年GCC邮件列表中发布的解释,GCC的最新版本处理中间计算中的过多精度.gcc -std=c99 -mno-sse2 -mfpmath=387据我所知,这个描述使程序编译完全可预测,到最后一点.如果它偶然没有,那就是一个错误而且它将被修复:约瑟夫S.迈尔斯在他的帖子中声明的意图是使其可预测.

是否记录了Clang如何处理超额精度(比如何时使用该选项-mno-sse2),以及在哪里？

(*)编辑:这是夸大其词.当允许将x87 FPU配置为使用53位有效数字时,这有点令人讨厌,但并不难以模拟binary64.

---

在下面的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)
…

Answer 1

这不能回答最初提出的问题,但如果您是处理类似问题的程序员,这个答案可能会对您有所帮助.

我真的不知道感知到的困难在哪里.提供严格的IEEE-754 binary64语义,同时限制为80387浮点数学,并保留80位长双倍计算,似乎遵循GCC-4.6.3和clang-3.0(基于LLVM)的明确指定的C99转换规则3.0).

编辑补充:然而,Pascal Cuoq是正确的:gcc-4.6.3或clang-llvm-3.0实际上都没有为'387浮点数学正确地执行这些规则.给定正确的编译器选项,规则正确应用于在编译时计算的表达式,但不适用于运行时表达式.有一些解决方法,在下面的休息后列出.

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

(我喜欢的一个优选示例是Kahan求和算法.使用C99和适当的转换(向维基百科示例代码添加强制转换),根本不需要任何技巧或额外的临时变量.无论编译器优化级别如何,实现都有效,包括-O3和-Ofast.)

C99明确指出(例如5.4.2.2),转换和赋值都删除了所有额外的范围和精度.这意味着您可以long double通过定义计算期间使用的临时变量来使用算术long double,同时将输入变量转换为该类型; 每当需要IEEE-754二进制64时,只需转换为a 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转换规则只是您必须学习的另一个工具.真的,你遇到的确切同样的问题,如果你混合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和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长期双倍(即ld == 0.1L),屈服true,而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都-mfpmath=387指定使用80位精度计算所有表达式.这导致例如

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

7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000

产生不正确的结果,因为二进制结果中间的那个字符串恰好位于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

但是,这只影响编译器优化的表达式,并且似乎根本不修复387处理.如果如使用clang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test与test.c被帕斯卡尔Cuoq的示例程序,你会得到每个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的编译器错误,也是一个容易复制的错误.(帕斯卡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>或open-code.例如,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需要扩展精度".幸运的是,您可以从GNU C库中获取'387实现,并在头文件中实现它们,或者libm如果您需要该math.h功能,则编写已知的工作; 事实上,我想我可以帮助那里.

Answer 2

为了记录,下面是我通过实验发现的.以下程序显示使用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  t.c && ./a.out
r1=0 r2=1 r3=0 r4=1 r5=1 r6=0 r7=1
$ clang -mno-sse2 -std=c99 -O2  t.c && ./a.out
r1=0 r2=1 r3=0 r4=1 r5=1 r6=1 r7=1

演员(double)区别于r5并r6在-O2具有没有影响-O0和变数r3和r4.结果r1与r5所有优化级别r6不同,而仅与r3at 不同-O2.