Ответ 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; на самом деле, я думаю, что мог бы помочь там.