Есть ли документ, описывающий, как Clang обрабатывает избыточную точность с плавающей точкой?
Почти невозможно (*) обеспечить строгую семантику IEEE 754 при разумных затратах, когда разрешено использовать только инструкции с плавающей запятой, 387. Это особенно трудно, когда хочется, чтобы FPU работал на полном 64-битном значении, и чтобы long double тип доступен для повышенной точности. Обычное "решение" заключается в выполнении промежуточных вычислений с единственно возможной точностью и преобразовании в более низкую точность в более или менее четко определенных случаях.
Последние версии GCC обрабатывают избыточную точность промежуточных вычислений в соответствии с интерпретацией, изложенной Джозефом С. Майерсом в публикации 2008 года в списке рассылки GCC. Это описание делает программу скомпилированной с gcc -std=c99 -mno-sse2 -mfpmath=387 насколько я понимаю, вполне предсказуемо, до последней черты. А если случайно это не так, то это ошибка, и она будет исправлена: заявленное Джозефом С. Майерсом намерение сделать его предсказуемым.
Задокументировано ли, как Clang обрабатывает избыточную точность (скажем, когда опция -mno-sse2 используется), а где?
(*) РЕДАКТИРОВАТЬ: это преувеличение. Немного раздражает, но не так сложно эмулировать двоичную 64-битную версию, когда разрешено конфигурировать x87 FPU для использования 53-битного значения.
После комментария 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)
…
2 ответа
Это не отвечает на изначально поставленный вопрос, но если вы программист, работающий с подобными проблемами, этот ответ может вам помочь.
Я действительно не вижу, где воспринимается трудность. Обеспечение строгой семантики IEEE-754 binary64 при ограничении математикой с плавающей запятой 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-битными значениями с плавающей запятой. Действительно, вы сталкиваетесь с точно такими же проблемами, если вы смешиваете поплавки binary32 и binary64 (float а также double на большинстве архитектур)!
Возможно, переписывание кода исследования Паскаля Куока, чтобы правильно применить правила приведения, делает это более понятным?
#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, поэтому они совершенно неинтересны для меня.
Паскаль Куок поднимает интересную проблему в цепочке комментариев ниже, которую я не сразу узнал.
При оценке выражения как GCC, так и лязг с -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 Будучи примером программы Паскаля Куока, вы получите правильный результат по правилам 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 функциональные возможности; на самом деле, я думаю, что смогу там помочь.
Для справки ниже приведено то, что я нашел экспериментально. Следующая программа показывает различные варианты поведения при компиляции с 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,