Различные результаты с плавающей запятой с включенной оптимизацией - ошибка компилятора?

Ответ 1

Процессоры Intel x86 используют внутреннюю 80-битную расширенную точность, тогда как double обычно 64-битная. Различные уровни оптимизации влияют на то, как часто значения с плавающей запятой из ЦП сохраняются в памяти и, таким образом, округляются от 80-битной точности до 64-битной точности.

Используйте параметр -ffloat-store gcc, чтобы получить одинаковые результаты с плавающей точкой с разными уровнями оптимизации.

В качестве альтернативы используйте тип long double, который обычно равен 80 бит на gcc, чтобы избежать округления от 80-битной до 64-разрядной точности.

man gcc говорит все:

   -ffloat-store
       Do not store floating point variables in registers, and inhibit
       other options that might change whether a floating point value is
       taken from a register or memory.

       This option prevents undesirable excess precision on machines such
       as the 68000 where the floating registers (of the 68881) keep more
       precision than a "double" is supposed to have.  Similarly for the
       x86 architecture.  For most programs, the excess precision does
       only good, but a few programs rely on the precise definition of
       IEEE floating point.  Use -ffloat-store for such programs, after
       modifying them to store all pertinent intermediate computations
       into variables.

Ответ 2

Выход должен быть: 4.5 4.6 Это будет результат, если у вас была бесконечная точность, или если вы работаете с устройством, использующим представление с плавающей запятой на основе десятичной системы, а не на основе двоичного кода. Но это не так. Большинство компьютеров используют двоичный стандарт IEEE с плавающей запятой.

Как уже отмечал Максим Егорушкин в своем ответе, часть проблемы заключается в том, что внутренне ваш компьютер использует 80-битное представление с плавающей запятой. Это лишь часть проблемы. Основой проблемы является то, что любое число формы n.nn5 не имеет точного двоичного плавающего представления. Эти угловые случаи всегда являются неточными числами.

Если вы действительно хотите, чтобы округление могло быть надежно повернуто вокруг этих угловых случаев, вам нужен алгоритм округления, который учитывает тот факт, что n.n5, n.nn5 или n.nnn5 и т.д. (но не n.5 ) всегда неточно. Найдите угловой регистр, который определяет, округляет ли или нет значение входного значения и возвращает округленное или округленное значение на основе сравнения с этим угловым случаем. И вам нужно позаботиться о том, чтобы оптимизирующий компилятор не помещал найденный угловой регистр в расширенный регистр точности.

См. Как Excel успешно объединяет плавающие числа, даже если они неточны? для такого алгоритма.

Или вы можете просто жить с тем, что угловые случаи иногда будут ошибочно округлены.

Ответ 3

Различные компиляторы имеют разные настройки оптимизации. Некоторые из этих более быстрых параметров оптимизации не поддерживают строгие правила с плавающей запятой в соответствии с IEEE 754. Visual Studio имеет определенную настройку /fp:strict, /fp:precise, /fp:fast, где /fp:fast нарушает стандарт на то, что можно сделать. Вы можете обнаружить, что этот флаг управляет оптимизацией в таких настройках. Вы также можете найти аналогичную настройку в GCC, которая изменяет поведение.

Если это так, то единственное, что отличается между компиляторами, это то, что GCC будет искать наиболее быстрое поведение с плавающей точкой по умолчанию при более высоких оптимизациях, тогда как Visual Studio не изменяет поведение с плавающей точкой с более высокими уровнями оптимизации. Таким образом, это может быть не обязательно реальная ошибка, но предполагаемое поведение варианта, о котором вы не знали, что вы включаете.

Ответ 4

Тем, кто не может воспроизвести ошибку: не раскомментируйте прокомментированные отладки stmts, они влияют на результат.

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

Дальнейший вопрос:

Мне было интересно, должен ли я всегда включать параметр -ffloat-store?

Чтобы быть легкомысленным, должна быть причина, по которой некоторые программисты не включаются -ffloat-store, иначе опция не будет существовать (также должна быть причина, по которой некоторые программисты делают включите -ffloat-store). Я бы не рекомендовал всегда включать его или всегда отключать. Включение этого параметра предотвращает некоторые оптимизации, но отключение этого режима позволяет получить поведение, которое вы получаете.

Но, как правило, существует некоторое несоответствие между двоичными числами с плавающей запятой (например, компьютер) и десятичными числами с плавающей запятой (что люди знакомы), и это несоответствие может привести к аналогичному поведению к тому, что вы получаете (чтобы быть ясным, поведение, которое вы получаете, не, вызванное этим несоответствием, но похожее поведение может быть). Дело в том, что, поскольку у вас уже есть неопределенность при работе с плавающей точкой, я не могу сказать, что -ffloat-store делает ее лучше или хуже.

Вместо этого вы можете захотеть взглянуть на другие решения на проблему, которую вы пытаетесь решить (к сожалению, Кениг не указывает на фактическую бумагу, и я не могу найти для нее явное "каноническое" место, поэтому мне придется отправить вас на Google).

Если вы не используете округление для вывода, я бы, вероятно, посмотрел на std::modf() (в cmath) и std::numeric_limits<double>::epsilon() (в limits). Думая о оригинальной функции round(), я считаю, что было бы проще заменить вызов std::floor(d + .5) вызовом этой функции:

// this still has the same problems as the original rounding function
int round_up(double d)
{
    // return value will be coerced to int, and truncated as expected
    // you can then assign the int to a double, if desired
    return d + 0.5;
}

Я думаю, это предполагает следующее улучшение:

// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
    double floor;
    d = std::modf(d, &floor);
    return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}

Простая заметка: std::numeric_limits<T>::epsilon() определяется как "наименьшее число, добавленное к 1, которое создает число, не равное 1." Обычно вам нужно использовать относительный epsilon (то есть, масштаб epsilon каким-то образом учитывает тот факт, что вы работаете с числами, отличными от "1" ). Сумма d, .5 и std::numeric_limits<double>::epsilon() должна быть около 1, поэтому группировка этого добавления означает, что std::numeric_limits<double>::epsilon() будет иметь нужный размер для того, что мы делаем. Если что-либо, std::numeric_limits<double>::epsilon() будет слишком большим (когда сумма всех трех меньше единицы) и может привести нас к округлению некоторых чисел, когда мы не должны.

В настоящее время вы должны рассмотреть std::nearbyint().

Ответ 5

Принятый ответ верен, если вы компилируете для цели x86, которая не включает SSE2. Все современные процессоры x86 поддерживают SSE2, поэтому, если вы можете воспользоваться этим, вам следует:

-mfpmath=sse -msse2 -ffp-contract=off

Давай сломаем это.

-mfpmath=sse -msse2. Это выполняет округление с использованием регистров SSE2, что намного быстрее, чем сохранение каждого промежуточного результата в памяти. Обратите внимание, что это уже значение по умолчанию для GCC для x86-64. Из вики GCC:

На более современных процессорах x86, которые поддерживают SSE2, указание опций компилятора -mfpmath=sse -msse2 гарантирует, что все операции с плавающей запятой и двойные операции выполняются в регистрах SSE и правильно округляются. Эти параметры не влияют на ABI и поэтому должны использоваться по возможности для предсказуемых числовых результатов.

-ffp-contract=off. Однако для точного совпадения недостаточно контроля округления. Команды FMA (слияния-умножения-сложения) могут изменить поведение округления по сравнению с его не слитыми аналогами, поэтому мы должны его отключить. Это значение по умолчанию для Clang, а не GCC. Как объясняется этим ответом:

FMA имеет только одно округление (оно эффективно сохраняет бесконечную точность для внутреннего результата временного умножения), в то время как ADD + MUL имеет два.

Отключая FMA, мы получаем результаты, которые точно совпадают при отладке и выпуске, за счет некоторой производительности (и точности). Мы все еще можем воспользоваться другими преимуществами производительности SSE и AVX.

Ответ 6

Я углубился в эту проблему, и я могу принести больше точности. Во-первых, точные представления 4.45 и 4.55 в соответствии с gcc на x84_64 следующие (с libquadmath для печати последней точности):

float 32:   4.44999980926513671875
double 64:  4.45000000000000017763568394002504646778106689453125
doublex 80: 4.449999999999999999826527652402319290558807551860809326171875
quad 128:   4.45000000000000000000000000000000015407439555097886824447823540679418548304813185723105561919510364532470703125

float 32:   4.55000019073486328125
double 64:  4.54999999999999982236431605997495353221893310546875
doublex 80: 4.550000000000000000173472347597680709441192448139190673828125
quad 128:   4.54999999999999999999999999999999984592560444902113175552176459320581451695186814276894438080489635467529296875

Как сказал Максим выше, проблема в 80-битном размере регистров FPU.

Но почему проблема никогда не происходит в Windows? на IA-32 FPU x87 был настроен на использование внутренней точности для мантиссы в 53 бита (эквивалентно общему размеру 64 бита: double). Для Linux и Mac OS была использована точность по умолчанию в 64 бита (эквивалентно общему размеру 80 бит: long double). Таким образом, проблема должна быть возможной или нет на этих разных платформах путем изменения управляющего слова FPU (при условии, что последовательность инструкций вызовет ошибку). О проблеме сообщили в gcc как об ошибке 323 (прочитайте хотя бы комментарий 92!).

Чтобы показать точность мантиссы в Windows, вы можете скомпилировать это в 32 бита с помощью VC++:

#include "stdafx.h"
#include <stdio.h>  
#include <float.h>  

int main(void)
{
    char t[] = { 64, 53, 24, -1 };
    unsigned int cw = _control87(0, 0);
    printf("mantissa is %d bits\n", t[(cw >> 16) & 3]);
}

и в Linux/Cygwin:

#include <stdio.h>

int main(int argc, char **argv)
{
    char t[] = { 24, -1, 53, 64 };
    unsigned int cw = 0;
    __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw));
    printf("mantissa is %d bits\n", t[(cw >> 8) & 3]);
}

Обратите внимание, что с помощью gcc вы можете установить точность FPU с помощью -mpc32/64/80, хотя в Cygwin это игнорируется. Но имейте в виду, что это изменит размер мантиссы, но не экспоненты, позволяя двери открыться другим видам другого поведения.

В архитектуре x86_64 SSE используется, как говорит tmandry, поэтому проблема не возникнет, если вы не принудительно используете старый FPU x87 для вычислений FP с -mfpmath=387 или если вы не скомпилируете в 32-битном режиме с -m32 (вы будете нужен мультибиблиотечный пакет). Я мог бы воспроизвести проблему в Linux с различными комбинациями флагов и версий gcc:

g++-5 -m32 floating.cpp -O1
g++-8 -mfpmath=387 floating.cpp -O1

Я попробовал несколько комбинаций на Windows или Cygwin с VC++/gcc/tcc, но ошибка так и не появилась. Я предполагаю, что последовательность сгенерированных инструкций не одинакова.

Наконец, обратите внимание, что экзотическим способом предотвращения этой проблемы с 4.45 или 4.55 было бы использование _Decimal32/64/128, но поддержка действительно скудна... Я потратил много времени, чтобы иметь возможность сделать printf с libdfp !

Ответ 7

Лично я столкнулся с той же проблемой, идя по другому пути - от gcc до VS. В большинстве случаев я считаю, что лучше избегать оптимизации. Единственный раз, когда стоит иметь дело с численными методами, связанными с большими массивами данных с плавающей запятой. Даже после разборки меня часто преследует выбор компиляторов. Очень часто просто проще использовать встроенные компиляторы или просто написать сборку самостоятельно.