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


приведенный ниже код работает на Visual Studio 2008 с оптимизацией и без нее. Но он работает только на g++ без оптимизации (O0).

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

вывод должен быть:

4.5
4.6

но g++ с оптимизацией (O1 -O3) будет выход:

4.5
4.5

если я добавлю volatile ключевое слово перед t, оно работает, так что может быть какая-то ошибка оптимизации?

тест на g++ 4.1.2 и 4.4.4.

вот результат на ideone: http://ideone.com/Rz937

и вариант, который я тестирую на g++, прост:

g++ -O2 round.cpp

более интересный результат, даже я включаю /fp:fast параметр в Visual Studio 2008, результат по-прежнему является правильным.

еще вопрос:

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

потому что версия g++, которую я тестировал,поставляется с CentOS/Red Hat Linux 5 и CentOS / Redhat 6.

я скомпилировал многие из моих программ под этими платформами, и я беспокоюсь, что это вызовет неожиданные ошибки внутри моих программ. Кажется немного сложным исследовать весь мой код C++ и используемые библиотеки, есть ли у них такие проблемы. Есть предложения?

кто-нибудь заинтересован в том, почему даже /fp:fast включен, Visual Studio 2008 все еще работает? Похоже, что Visual Studio 2008-это более надежный в этой проблеме, чем g++?

7 98

7 ответов:

процессоры 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.

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

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

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

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

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

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

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

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

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

еще вопрос:

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

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

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

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


если вы не округляете для целей вывода, я бы, вероятно, посмотрел на std::modf() (in cmath) и std::numeric_limits<double>::epsilon() (in 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."Обычно вам нужно использовать относительный Эпсилон (т. е. масштабировать Эпсилон каким-то образом, чтобы учесть тот факт, что вы работаете цифры, кроме "1"). Сумма d,.5 и std::numeric_limits<double>::epsilon() должно быть около 1, поэтому группировка, что добавление означает, что std::numeric_limits<double>::epsilon() будет правильный размер для того, что мы делаем. Во всяком случае,std::numeric_limits<double>::epsilon() будет слишком большим (когда сумма всех трех меньше единицы) и может заставить нас округлить некоторые числа, когда мы не должны.


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

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

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

давайте разберем это.

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

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

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

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

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

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

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

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, так что проблема не произойдет, если вы не заставите старый x87 FPU для вычисления FP с помощью -mfpmath=387, или если вы не компилируете в 32-битном режиме с -m32 (вам понадобится пакет multilib). Я мог бы воспроизвести проблему на 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 !