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

Question

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

根据 Joseph S. Myers 在 2008 post to the GCC mailing list 中提出的解释，最新版本的 GCC 处理中间计算中的超额精度。这个描述使得用 gcc -std=c99 -mno-sse2 -mfpmath=387 编译的程序完全可以预测，据我所知。如果碰巧没有，它就是一个错误，它将被修复:Joseph S. Myers 在他的帖子中声明的意图是使其可预测。

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

(*) 编辑:这是夸大其词。当允许将 x87 FPU 配置为使用 53 位有效数时，模拟 binary64 为 slightly annoying but not that difficult。

在下面 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 位长 double 计算，似乎遵循 GCC-4.6.3 和 clang-3.0(基于 LLVM 3.0)。

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

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

(我喜欢的一个首选示例是 Kahan 求和算法。使用 C99 和适当的转换(将转换添加到例如 Wikipedia 示例代码)，根本不需要技巧或额外的临时变量。无论编译器优化级别如何，实现都可以工作，包括在-O3 和 -Ofast .)

C99 明确声明(例如 5.4.2.2)强制转换和赋值都删除了所有额外的范围和精度。这意味着您可以通过将计算期间使用的临时变量定义为 0x2518122231343141 来使用 long double 算术，并将输入变量转换为该类型；每当需要 IEEE-754 binary64 时，只需将其转换为 long double 。

在 '387 上，强制转换在上述两个编译器上生成赋值和加载；这确实将 80 位值正确舍入为 IEEE-754 binary64。这个成本在我看来是非常合理的。所需的确切时间取决于体系结构和周围的代码；通常它是并且可以与其他代码交织以将成本降低到可以忽略的水平。当 MMX、SSE 或 AVX 可用时，它们的寄存器与 80 位 80387 寄存器是分开的，并且通常通过将值移动到 MMX/SSE/AVX 寄存器来完成转换。

(我更喜欢生产代码使用特定的浮点类型，比如 double 等，用于临时变量，以便它可以定义为 tempdouble 或 double 取决于所需的架构和速度/精度权衡)。

简而言之:

Don't assume (expression) is of double precision just because all the variables and literal constants are. Write it as (double)(expression) if you want the result at double precision.

这也适用于复合表达式，并且有时可能会导致具有多个转换级别的笨拙的表达式。

如果您希望以 80 位精度计算 long double 和 expr1，但还需要先将每个四舍五入到 64 位的乘积，请使用

long double  expr1;
long double  expr2;
double       product = (double)(expr1) * (double)(expr2);

注意， expr2 计算为两个 64 位值的乘积；不以 80 位精度计算，然后向下舍入。以 80 位精度计算乘积，然后四舍五入，将是

double       other = expr1 * expr2;

或者，添加描述性类型转换，告诉你到底发生了什么，

double       other = (double)((long double)(expr1) * (long double)(expr2));

很明显， product 和 product 经常不同。

如果您确实使用混合的 32 位/64 位/80 位/128 位浮点值，C99 转换规则只是您必须学会使用的另一个工具。真的，如果您混合使用 binary32 和 binary64 浮点数(在大多数架构上为 other 和 float)，您会遇到完全相同的问题!

也许重写 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 将 double 的右侧提升为 long double first(即到 ld == 0.1 )，产生 ld == 0.1L ，而对于 SSE/AVX，则为 0x214138131313

对于纯 '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

各种优化标志组合为 true ，包括 long double 、 ... 、 -fomit-frame-pointer 、 0x251812231318141 、 0x2518143 13 13 13 13 13 13 13 14 341 41 2

使用任何其他标志或 C99 编译器应该会导致相同的结果，除了 -O0 大小(以及当前 GCC 版本的 -O1)。如果您遇到任何差异，我将非常感激听到它们；我可能需要警告我的用户使用此类编译器(编译器版本)。请注意，Microsoft 不支持 C99，因此我对它们完全不感兴趣。

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

在计算表达式时，GCC 和带有 -O2 的 clang 指定所有表达式都使用 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

仅用 -O3 获得的结果是

42695510550671098263984292201741942784 = 0x1.00f6d98d0a43p+125 = 100000000111101101101100110001101000010100100001100000000000000000000000000000000000000000000000000000000000000000000000000000

理论上，您应该能够通过使用选项告诉 gcc 和 clang 强制执行正确的 C99 舍入规则

-std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard

但是，这只影响编译器优化的表达式，似乎根本无法修复 387 处理。如果您使用例如 -Os 与 long double 为 Pascal Cuoq's example program ，您将根据 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 指出 ld == 1.0 意味着对 double 参数的操作总是以扩展精度完成，但我看不出有任何合理的理由——除了必须重写 '387 上的 -mfpmath=387 的一部分——在 C99 中这样做并考虑IEEE-754 规则可由硬件实现!毕竟，通过修改 '387 控制字以匹配表达式的精度，编译器很容易实现正确的操作。而且，考虑到应该强制这种行为的编译器选项 - - -std=c99 -m32 -mno-sse -mfpmath=387 -- 如果实际上需要 clang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test 行为没有意义，也不存在向后兼容性问题。)重要的是要注意，给定正确的编译器标志，在编译时计算的表达式确实得到正确计算，并且只有表达式在运行时评估得到不正确的结果。

最简单的解决方法，也是可移植的方法，是使用 test.c(来自 FLT_EVAL_METHOD=2)将所有结果四舍五入到零。

在某些情况下，向零舍入可能有助于可预测性和病理情况。特别是，对于像 libm 这样的区间，向零舍入意味着永远不会通过舍入达到上限；重要的是，如果您评估例如分段样条。

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

您可以使用来自 -std=c99 -ffloat-store -fexcess-precision=standard 的 FLT_EVAL_METHOD=2 或对其进行开放编码。例如， fesetround(FE_TOWARDZERO) :

#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

换句话说，您可以通过使用 fenv.h 设置精度和舍入模式来解决编译器问题。

缺点是某些数学库( x = [0,1) 、 __FPU_SETCW() )可能会假设中间结果始终以 80 位精度计算。至少 x86_64 上的 GNU C 库 #include <fpu_control.h> 有注释“libm 需要扩展精度”。幸运的是，您可以从例如 '387 实现。 GNU C 库，如果需要 precision.c 功能，则在头文件中实现它们或编写已知的工作 fp387() ；事实上，我想我或许能帮上忙。