Visual C ++ math.h ошибка

Question

Я отлаживал свой проект и не смог найти ошибку. Наконец я нашел это. Посмотри на код. Вы думаете, что все в порядке, и результат будет "OK! OK! OK!", Не так ли? Теперь скомпилируйте его с помощью VC (я пробовал vs2005 и vs2008).

#include <math.h>
#include <stdio.h>


int main () {
    for ( double x = 90100.0; x<90120.0; x+=1 )
    {
        if ( cos(x) == cos(x) )
            printf ("x==%f  OK!\n", x);
        else
            printf ("x==%f  FAIL!\n", x);
    }

    getchar();
    return 0; 
}

Магическая двойная константа 90112.0. Когда x <90112.0, все в порядке, когда x> 90112.0 - Нет! Вы можете изменить cos на грех.

Есть идеи? Не забывайте, что грех и cos периодичны.

Answer 1

Может быть так: http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18

Я знаю, что это трудно принять, но арифметика с плавающей запятой просто не работает, как большинство людей ожидают. Хуже того, некоторые различия зависят от деталей аппаратного обеспечения вашего компьютера с плавающей запятой и / или настроек оптимизации, которые вы используете в своем конкретном компиляторе. Тебе это может не понравиться, но так оно и есть. Единственный способ «получить это» - отложить свои предположения о том, как вещи должны вести себя и принимать вещи так, как они на самом деле делают ведут себя ...

(с акцентом на слово «часто»; поведение зависит от вашего оборудования, компилятора и т. Д.): Вычисления и сравнения с плавающей запятой часто выполняются специальным оборудованием, которое часто содержит специальные регистры, и эти регистры часто имеют больше битов чем double. Это означает, что промежуточные вычисления с плавающей запятой часто имеют больше битов, чем sizeof(double), и когда значение с плавающей запятой записывается в ОЗУ, оно часто усекается, часто теряя некоторые биты точности ...

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

Answer 2

Как уже отмечали другие, математическая библиотека VS выполняет вычисления на FPU x87 и генерирует 80-битные результаты, даже если тип имеет тип double.

Таким образом:

1. вызывается cos () и возвращается с cos (x) на вершине стека x87 в виде 80-битного числа с плавающей запятой
2. cos (x) извлекается из стека x87 и сохраняется в памяти как double; это приводит к его округлению до 64-битного числа с плавающей запятой, что меняет свое значение
3. вызывается cos () и возвращается с cos (x) на вершине стека x87 в виде 80-битного числа с плавающей запятой
4. округленное значение загружается в стек x87 из памяти
5. округленные и необоснованные значения cos (x) сравниваются неравно.

Многие математические библиотеки и компиляторы защищают вас от этого либо путем выполнения вычислений в 64-битных числах с плавающей запятой в регистрах SSE, когда они доступны, либо путем принудительного хранения и округления значений перед сравнением, либо путем сохранения и перезагрузки конечного результата фактическое вычисление cos (). Комбинация компилятор / библиотека, с которой вы работаете, не так проста.

Answer 3

Процедура cos (x) == cos (x), сгенерированная в режиме разблокировки:

00DB101A  call        _CIcos (0DB1870h) 
00DB101F  fld         st(0) 
00DB1021  fucompp 

Значение вычисляется один раз, затем клонируется, затем сравнивается с самим собой - результат будет в порядке

То же самое в режиме отладки:

00A51405  sub         esp,8 
00A51408  fld         qword ptr [x] 
00A5140B  fstp        qword ptr [esp] 
00A5140E  call        @ILT+270(_cos) (0A51113h) 
00A51413  fld         qword ptr [x] 
00A51416  fstp        qword ptr [esp] 
00A51419  fstp        qword ptr [ebp-0D8h] 
00A5141F  call        @ILT+270(_cos) (0A51113h) 
00A51424  add         esp,8 
00A51427  fld         qword ptr [ebp-0D8h] 
00A5142D  fucompp          

Теперь происходят странные вещи.
1. X загружен в fstack (X, 0)
2. Х хранится в обычном стеке (усечение)
3. Косинус вычисляется, результат на стеке с плавающей точкой
4. Х загружается снова
5. X хранится в обычном стеке (усечение, так как на данный момент мы «симметричны»)
6. Результат 1-го косинуса, который был в стеке, сохраняется в памяти, теперь происходит еще одно усечение для 1-го значения
7. Косинус вычисляется, 2-й результат, если он находится в стеке float, но это значение было усечено только один раз
8. 1-е значение загружается в fstack, но это значение было усечено дважды (один раз до вычисления косинуса, один раз после)
9. Эти два значения сравниваются - мы получаем ошибки округления.

Answer 4

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

Регистры с плавающей запятой могут иметь размер, отличающийся от значений памяти (в современных машинах Intel регистры FPU имеют значение 80 бит против 64 бит). Если компилятор генерирует код, который вычисляет первый косинус, затем сохраняет значение в памяти, вычисляет второй косинус и сравнивает значение в памяти с значением в регистре, тогда значения могут отличаться (из-за проблем округления от 80 до 64 бит) ,

Значения с плавающей точкой немного хитры. Google для сравнения с плавающей запятой.

Answer 5

Компилятор мог сгенерировать код, который заканчивает тем, что сравнивал 64-битное двойное значение с 80-битный внутренний регистр с плавающей запятой. Проверка значений с плавающей запятой на равенство склонен к такого рода ошибкам - вам почти всегда лучше делать «нечеткое» сравнение, например (fabs (val1 - val2)

Answer 6

Увеличение и проверка значения с плавающей запятой как переменной управления цикла, как правило, является действительно плохой идеей. Создайте отдельный int LCV только для зацикливания, если вам нужно.

В этом случае проще:

for ( int i = 90100; i<90120; i+=1 )    {
    if ( cos(i) == cos(i) )
        printf ("i==%d  OK!\n", i);
    else
        printf ("i==%d  FAIL!\n", i);
}

Answer 7

Как обойти проблему? Изменить , если блок:

if ( (float)cos(x) == (float)cos(x) )