Есть ли документ, описывающий, как Clang обрабатывает избыточную точность с плавающей запятой?

Ответ 1

Это не отвечает первоначально поставленному вопросу, но если вы программист, работающий с подобными проблемами, этот ответ может вам помочь.

Я действительно не вижу, где воспринимается трудность. Предоставление строгой семантики binary64 IEEE-754 при ограничении математикой с плавающей запятой 80387 и сохранение 80-битного двойного вычисления, по-видимому, соответствует хорошо определенным правилам каста C99 как с GCC-4.6.3, так и с clang-3.0 (на основе LLVM 3.0).

Отредактировано для добавления: Тем не менее, Pascal Cuoq верен: ни gcc-4.6.3, ни clang-llvm-3.0 не применяют эти правила правильно для математики с плавающей запятой 387. Учитывая правильные параметры компилятора, правила корректно применяются к выражениям, оцениваемым во время компиляции, но не для выражений во время выполнения. Возможны обходные пути, перечисленные после перерыва ниже.

Я использую код моделирования молекулярной динамики и очень хорошо знаком с требованиями к повторяемости/предсказуемости, а также с желанием сохранить максимальную точность, когда это возможно, поэтому я утверждаю, что знаю, о чем я говорю. Этот ответ должен показать, что инструменты существуют и просты в использовании; проблемы возникают из-за того, что они не знают или не используют эти инструменты.

(Предпочтительный пример, который мне нравится, это алгоритм суммирования Кахана. С C99 и правильной литой (добавление отливок, например, в код примера Википедии) вообще не нужны никакие трюки или дополнительные временные переменные. Реализация работает независимо от уровня оптимизации компилятора, включая -O3 и -Ofast.)

C99 явно заявляет (например, 5.4.2.2), что литье и присвоение удаляют весь дополнительный диапазон и точность. Это означает, что вы можете использовать арифметику long double, определяя ваши временные переменные, используемые при вычислении, как long double, также применяя ваши входные переменные к этому типу; всякий раз, когда требуется бинарный бит IEEE-754, просто добавьте к double.

В '387, литье генерирует назначение и нагрузку на обоих вышеперечисленных компиляторах; это корректно соответствует 80-битовому значению для бинарного кода IEEE-754. На мой взгляд, это очень разумно. Точное время зависит от архитектуры и окружающего кода; обычно это и может чередоваться с другим кодом, чтобы снизить стоимость до уровня пренебрежимости. Когда доступны 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 часто различаются.

Правила кастинга C99 - это еще один инструмент, который вы должны научиться владеть, если вы работаете со смешанными 32-битными/64-битными/80-битными/128-битными значениями с плавающей запятой. Действительно, вы сталкиваетесь с теми же проблемами, если вы смешиваете бинарные32 и бинарные64 float (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 size (и ld == 1.0 для текущих версий GCC). Если у вас возникнут какие-либо разногласия, я буду очень благодарен вам за это; Возможно, мне придется предупредить моих пользователей о таких компиляторах (версии компилятора). Обратите внимание, что Microsoft не поддерживает C99, поэтому они мне совершенно неинтересны.

Паскаль Куок вызывает интересную проблему в цепочке комментариев ниже, о которой я не сразу узнал.

При оценке выражения как 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 пример программы Pascal 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, и легко воспроизведенную. (Паскаль Куок указывает, что FLT_EVAL_METHOD=2 означает, что операции с аргументами с двойной точностью всегда выполняются с расширенной точностью, но я не вижу разумной причины - кроме необходимости переписывать части libm на '387 - делать это на C99 и учитывая, что правила IEEE-754 достижимы аппаратными средствами! В конце концов, правильная операция легко достижима компилятором путем изменения управляющего слова "387" в соответствии с точностью выражения. И, учитывая параметры компилятора, которые должны заставить это поведение - -std=c99 -ffloat-store -fexcess-precision=standard - не имеет смысла, если поведение FLT_EVAL_METHOD=2 действительно желательно, также нет проблем с обратной совместимостью.) Важно отметить, что, учитывая правильные флаги компилятора, выражения, оцененные во время компиляции, получить оценку правильно и что только выражения, оцененные во время выполнения, получают неверные результаты.

Простейшим обходным решением и переносным является использование fesetround(FE_TOWARDZERO) (от fenv.h) для округления всех результатов до нуля.

В некоторых случаях округление к нулю может помочь в предсказуемости и патологических случаях. В частности, для интервалов, подобных x = [0,1), округление к нулю означает, что верхний предел никогда не достигается округлением; важно, если вы оцениваете, например. кусочно сплайны.

Для других режимов округления вам необходимо напрямую управлять оборудованием 387.

Вы можете использовать либо __FPU_SETCW() от #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 бит. По крайней мере, библиотека GNU C fpu_control.h на x86_64 имеет комментарий "libm требует расширенной точности". К счастью, вы можете использовать "387 реализаций", например, Библиотеку GNU C и реализовать их в файле заголовка или записать известный для работы libm, если вам нужна функциональность math.h; на самом деле, я думаю, что мог бы помочь там.

Ответ 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 отличается от r3 на -O2.