Причудливое поведение с плавающей запятой с или без дополнительных переменных, почему?

Когда я запускаю следующий код в VC ++ 2013 (32-разрядный, без оптимизации):

#include <cmath>
#include <iostream>
#include <limits>

double mulpow10(double const value, int const pow10)
{
    static double const table[] =
    {
        1E+000, 1E+001, 1E+002, 1E+003, 1E+004, 1E+005, 1E+006, 1E+007,
        1E+008, 1E+009, 1E+010, 1E+011, 1E+012, 1E+013, 1E+014, 1E+015,
        1E+016, 1E+017, 1E+018, 1E+019,
    };
    return pow10 < 0 ? value / table[-pow10] : value * table[+pow10];
}

int main(void)
{
    double d = 9710908999.008999;
    int j_max = std::numeric_limits<double>::max_digits10;
    while (j_max > 0 && (
        static_cast<double>(
            static_cast<unsigned long long>(
                mulpow10(d, j_max))) != mulpow10(d, j_max)))
    {
        --j_max;
    }
    double x = std::floor(d * 1.0E9);
    unsigned long long y1 = x;
    unsigned long long y2 = std::floor(d * 1.0E9);
    std::cout
        << "x == " << x << std::endl
        << "y1 == " << y1 << std::endl
        << "y2 == " << y2 << std::endl;
}

я получил

x  == 9.7109089990089994e+018
y1 == 9710908999008999424
y2 == 9223372036854775808

в отладчике.

Я ошеломлен. Может кто-нибудь объяснить мне, как чертовски y1 и y2 имеют разные значения?

Обновлять:

Это происходит только в /Arch:SSE2 или /Arch:AVX, но не в /Arch:IA32 или /Arch:SSE.


c++ double floating-point visual-c++ unsigned-long-long-int
person user541686    schedule 31.01.2014    source источник
comment
Разница в том, что в одном случае это 64-битное число с плавающей запятой, а в другом - 80-битное число с плавающей запятой (именно так представлено число с плавающей запятой в математическом сопроцессоре). Это должно дать вам хорошее начало.   -  person sharptooth    schedule 31.01.2014
comment
@sharptooth: Извините, что? Разве здесь все не double, то есть 64-битное?   -  person user541686    schedule 31.01.2014
comment
@Mehrdad Арифметика с плавающей запятой может выполняться с большей точностью, чем указано ее типом, и ситуации, в которых это происходит, не нуждаются в документировании.   -  person    schedule 31.01.2014
comment
@ Mehrdad: Больше я не знаю подробностей, просто помню, что слышал что-то в этом роде. Я бы написал ответ, если бы у меня было лучшее объяснение.   -  person sharptooth    schedule 31.01.2014
comment
@hvd: Эээ, для меня это не имеет смысла - значение x, которое я вижу, является правильным, независимо от того, выполнялась ли его арифметика с 64-битной или 80-битной точностью. Поскольку это 64-битная переменная, это означает, что для ее представления было достаточно 64 бит. Это также означает, что преобразование его в unsigned long long должно дать правильный ответ, 64 или 80 бит. Кого волнует, является ли x 64-битным или 80-битным, если оно имеет правильное значение? Как это на что-нибудь влияет?   -  person user541686    schedule 31.01.2014
comment
@Mehrdad По какой логике это правильно? Это не равно математическому значению std::floor(9710908999.0089989 * 1.0E9). У вас есть точка зрения, что значение y2 слишком сильно отличается, чтобы быть вызванным этим, но это нормально, если оно не точно равно y1.   -  person    schedule 31.01.2014
comment
@hvd: Да, я знаю, что с плавающей запятой неточно, вы полностью упускаете суть, придирчиво. Под правильным я просто подразумеваю, что это имеет смысл, то есть в пределах ошибки округления. 9223372036854775808 неверно или в пределах ошибки округления, это признак некоторого переполнения.   -  person user541686    schedule 31.01.2014
comment
Как бы то ни было, я не воспроизводю в Linux с помощью gcc, который, как я знаю, имеет проблемы с повышенной точностью.   -  person AProgrammer    schedule 31.01.2014
comment
@Mehrdad Согласен, это не похоже на разумный результат.   -  person    schedule 31.01.2014
comment
Другие компиляторы кажутся более последовательными: codepad.org/fuoERz3k дает y1 = y2 = 9710908999008999424.   -  person Paul R    schedule 31.01.2014
comment
@AProgrammer: Смотрите мое обновление. Он отлично работает и в VC ++, если я просто не использую SSE2. Так что я действительно сбит с толку ... это ошибка кодогенератора? Это происходит в режиме отладки ...   -  person user541686    schedule 31.01.2014
comment
is it a codegen bug? Я бы не стал искать другого объяснения. Ошибки с плавающей запятой обычно возникают в компиляторах Microsoft. На одной из моих работ нам пришлось перевести несколько модулей на icc из-за них.   -  person n. 1.8e9-where's-my-share m.    schedule 31.01.2014
comment
Учитывая небольшой размер примера, я бы просто посмотрел на сгенерированный код, чтобы узнать, что происходит.   -  person JasonD    schedule 31.01.2014
comment
@Mehrdad Какая у вас версия компилятора VS? После тестирования с VS2013 он показывает согласованные результаты с каждой настройкой дуги.   -  person galop1n    schedule 31.01.2014
comment
Связываете неправильную стандартную библиотеку? (Я не пользователь VS, даже не знаю, есть ли у них разные стандартные библиотеки для параметров SSE / SSE2)   -  person AProgrammer    schedule 31.01.2014
comment
@ galop1n: Странно. У меня версия 18.00.21005.1   -  person user541686    schedule 31.01.2014
comment
@AProgrammer: Я так не думаю? Это должен быть стандарт.   -  person user541686    schedule 31.01.2014
comment
@ n.m .: Интересно, я немного не решаюсь назвать это ошибкой, мне интересно, может ли это быть разрешено C ++. Я выложу разборку через несколько минут.   -  person user541686    schedule 31.01.2014
comment
@Mehrdad Я не смог воспроизвести его с помощью VC ++ (2012), независимо от аргументов /Arch.   -  person James Kanze    schedule 31.01.2014
comment
Я просто понял, что это, похоже, зависит от кода вокруг него (я тестировал его в середине другой функции) ... позвольте мне посмотреть, могу ли я привести отдельный пример.   -  person user541686    schedule 31.01.2014
comment
Кто-нибудь еще указал, что 9223372036854775808 == 2 ^ 63? Мне кажется, это остатки жучка.   -  person laune    schedule 31.01.2014
comment
@JamesKanze: Я обновил его фрагментом кода, чтобы воспроизвести его. Вы можете воспроизвести это?   -  person user541686    schedule 31.01.2014
comment
@ galop1n: Я обновил вопрос и дал образец кода, который воспроизводит ошибку для меня.   -  person user541686    schedule 31.01.2014
comment
@ Mehrdad Я до сих пор не могу воспроизвести это. Я добавил вывод в ваш код, чтобы видеть значения. Но, возможно, это зависит от вариантов оптимизации.   -  person James Kanze    schedule 31.01.2014
comment
Я получаю тот же результат, что и вы. Я вижу вызов __ftol2 (преобразовать в long long со знаком), который соответствует возможной проблеме, предложенной в ответах, но не имеет для меня никакого смысла.   -  person    schedule 31.01.2014
comment
@JamesKanze: На самом деле у меня не включена какая-либо оптимизация. Флаги компилятора "%VS120COMNTOOLS%..\..\VC\bin\cl.exe" /nologo /I "%VS120COMNTOOLS%..\..\VC\include" Temp.cpp /link /LibPath:"%VS120COMNTOOLS%..\..\VC\lib" /LibPath:"%ProgramFiles(x86)%\Microsoft SDKs\Windows\v7.1A\Lib" запускают это нормально. Вы печатаете через std::cout или через printf? Почему-то кажется, что printf не работает, он печатает нули для меня.   -  person user541686    schedule 31.01.2014
comment
@hvd: Спасибо, рад узнать, что вы смогли воспроизвести это. Странная вещь заключается в следующем: если вы избавитесь от условия static_cast<double>(static_cast<unsigned long long>(mulpow10(d, j_max))) != mulpow10(d, j_max) в цикле, вы больше не получите этот результат, даже если он все еще используется __ftol2!   -  person user541686    schedule 31.01.2014
comment
Хм. __ftol2 фактически используется как для преобразования в long long со знаком, так и в long long без знака, но при преобразовании в unsigned long long его окружает довольно некоторый дополнительный код для обработки значений, выходящих за пределы допустимого диапазона.   -  person    schedule 31.01.2014
comment
@hvd: Похоже, позвонив _fpreset() до того, как y2 = std::floor(d * 1.0E9) исправит это, я не понимаю, почему.   -  person user541686    schedule 31.01.2014
comment
@Mehrdad Я смотрел сгенерированный код, и в нерабочей версии правильное значение хранится в памяти, другое значение загружается из памяти, а затем только что сохраненное (!) Значение снова загружается из памяти, но это не дает первоначального значения.   -  person    schedule 31.01.2014
comment
@hvd: смеется. Не знаю, что сказать ...   -  person user541686    schedule 31.01.2014
comment
@Mehrdad Похоже, произошло переполнение стека FPU. Некоторая часть цикла while оставляет значения в стеке FPU, и загрузка только что сохраненного значения позже не выполняется, потому что для этого значения больше нет места.   -  person    schedule 31.01.2014
comment
@hvd: Хм, что могло быть причиной этого? К сожалению, я не знаю, как работает FPU, зачем ему стек и т. Д.   -  person user541686    schedule 31.01.2014
comment
@Mehrdad FPU имеет восемь регистров, но использует их в виде стека, поэтому загрузка значения с плавающей запятой в регистр всегда загружает его в регистр ST0. Другие регистры перемещаются: ST1 получает старое значение ST0, ST2 получает старое значение ST1 и т. Д. Существуют также команды выталкивания для освобождения стека. Операции push и pop должны быть сбалансированы, иначе что-то сломается. :)   -  person    schedule 01.02.2014
comment
@hvd: Интересно, спасибо за информацию. Я немного сбит с толку, почему этого не происходит, когда я указываю / Arch: IA32?   -  person user541686    schedule 01.02.2014
comment
Эта проблема возникает в реализации вспомогательной функции преобразования в SSE2.   -  person    schedule 01.02.2014
comment
@hvd: Значит, в реализации SSE2 есть ошибка FPU x87?   -  person user541686    schedule 01.02.2014
comment
Да, он использует как FPU инструкции, так и инструкции SSE2 :)   -  person    schedule 01.02.2014

Ответы (4)


arrow_upward
5
arrow_downward

Вы конвертируете значения, выходящие за пределы double, в unsigned long long. Это недопустимо в стандартном C ++, и Visual C ++, похоже, действительно плохо обращается с этим в режиме SSE2: он оставляет число в стеке FPU, в конечном итоге переполняя его и приводя к сбоям в последующем коде, использующем FPU. действительно интересные способы.

Уменьшенный образец

double d = 1E20;
unsigned long long ull[] = { d, d, d, d, d, d, d, d };
if (floor(d) != floor(d)) abort();

Это прерывается, если ull имеет восемь или более элементов, но проходит, если у него до семи.

Решение состоит в том, чтобы не преобразовывать значения с плавающей запятой в целочисленный тип, если вы не знаете, что значение находится в диапазоне.

4.9 Преобразования с плавающей запятой [conv.fpint]

Prvalue типа с плавающей запятой можно преобразовать в prvalue целочисленного типа. Преобразование обрезается; то есть дробная часть отбрасывается. Поведение не определено, если усеченное значение не может быть представлено в целевом типе. [Примечание. Если тип назначения - bool, см. 4.12. - конец примечания]

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

Как бы то ни было, это не кажется преднамеренным, поэтому, хотя стандарт разрешает такое поведение, все же, возможно, стоит сообщить об этом как об ошибке.

person Community    schedule 31.01.2014
comment
+1 вау, я не знал, что это неопределенное поведение. Спасибо за анализ! - person user541686; 01.02.2014
comment
Кстати, как мне проверить, находится ли значение в пределах допустимого диапазона? Я не могу преобразовать std::numeric_limits<unsigned long long>::max() в double, так как он потеряет точность ... - person user541686; 01.02.2014
comment
@Mehrdad Я думаю, вы могли бы использовать < max() для общего случая и добавить специальную обработку для крайнего случая, но я не совсем уверен. - person ; 01.02.2014
comment
@Mehrdad Если подумать еще раз, вам нужно знать, является ли double(std::numeric_limits<unsigned long long>::max()) точным, округленным в меньшую или большую сторону. Я не уверен, как это сделать проще всего, но в крайнем случае преобразуйте его в строку. Если это точно, то преобразование безопасно, если x - double(std::numeric_limits<unsigned long long>::max()) < 1. Если округлить, то преобразование безопасно, если x < double(std::numeric_limits<unsigned long long>::max()). Если оно округлено в меньшую сторону, то преобразование безопасно, если x <= std::numeric_limits<unsigned long long>::max(). - person ; 01.02.2014
comment
@Mehrdad Я надеюсь, что есть способ написать это простым способом, который работает во всех трех случаях, но я не могу ничего придумать. (Я игнорирую отрицательные значения, как вы, наверное, заметили.) - person ; 01.02.2014
comment
Однако он переходит в unsigned long long, а не в signed long long. Результат находится в пределах диапазона unsigned long long. - person tmyklebu; 01.02.2014
comment
@tmyklebu Преобразование, завершившееся ошибкой, находится в пределах допустимого диапазона, но более раннее преобразование - нет, и проблема заключается в том, что более раннее преобразование является проблемой. - person ; 01.02.2014

arrow_upward
5
arrow_downward

9223372036854775808 равно 0x8000000000000000; то есть он равен INT64_MIN преобразованию к uint64_t.

Похоже, ваш компилятор преобразует возвращаемое значение floor в long long, а затем преобразует этот результат в unsigned long long.

Обратите внимание, что переполнение при преобразовании числа с плавающей запятой в целое является обычным делом, чтобы получить наименее представимое значение (например, cvttsd2siq на x86-64):

Если преобразование неточное, возвращается усеченный результат. Если преобразованный результат превышает максимальное целое число двойного слова со знаком, возникает недопустимое исключение с плавающей запятой, и если это исключение замаскировано, возвращается неопределенное целочисленное значение (80000000H).

(это из документации по двойному слову, но поведение четверного слова такое же.)

person ecatmur    schedule 31.01.2014
comment
Почему компилятор должен преобразовывать std::floor(9710908999.0089989 * 1.0E9) в long long, а затем в unsigned long long, но не преобразовывать x в long long, а затем в unsigned long long? - person Eric Postpischil; 31.01.2014
comment
Интересно, что я не получаю никаких исключений FP, даже когда использую _1 _..., но хорошее замечание о значении. - person user541686; 31.01.2014

arrow_upward
2
arrow_downward

Гипотеза: это ошибка. Компилятор правильно преобразует double в unsigned long long, но неправильно преобразует числа с плавающей запятой повышенной точности (возможно, long double) в unsigned long long. Подробности:

double              x = std::floor(9710908999.0089989 * 1.0E9);

Это вычисляет значение справа и сохраняет его в x. Значение в правой части может быть вычислено с повышенной точностью, но оно, как того требуют правила C ++, преобразуется в double при сохранении в x. Точное математическое значение будет 9710908999008998870, но округление его до формата double дает 9710908999008999424.

unsigned long long y1 = x;

Это преобразует значение double в x в unsigned long long, получая ожидаемое 9710908999008999424.

unsigned long long y2 = std::floor(9710908999.0089989 * 1.0E9);

Это вычисляет значение в правой части с использованием повышенной точности, получая 9710908999008998870. Когда значение повышенной точности преобразуется в unsigned long long, возникает ошибка, приводящая к 2 63 (9223372036854775808). Это значение, вероятно, является значением ошибки «вне допустимого диапазона», созданным инструкцией, преобразующей формат повышенной точности в 64-битное целое число. Компилятор использовал неправильную последовательность инструкций для преобразования своего формата повышенной точности в unsigned long long.

person Eric Postpischil    schedule 31.01.2014
comment
Я не думаю, что этот ответ правильный, потому что он отлично работает на / Arch: IA32 ... не должен ли он сломаться, если у него были проблемы с повышенной точностью? - person user541686; 31.01.2014
comment
+1 logtwo (9710908999.0089989 * 1.0E9) = 63.0743 явно вписывается в беззнаковый длинный - person ; 31.01.2014
comment
@Mehrdad Разные платформы / разные последовательности инструкций - person ; 31.01.2014
comment
@Mehrdad: Код компилятора для генерации инструкций для IA-32 должен несколько отличаться от кода для генерации инструкций для Intel 64. Код может быть написан разными людьми с разницей в годы. Конечно, у одного может быть ошибка, а у другого нет. - person Eric Postpischil; 31.01.2014
comment
@EricPostpischil: Что такое Intel 64? Я не тестирую на Intel 64, это весь 32-битный код. Я тестирую только / Arch: IA32 и / Arch: SSE2. Но если вы думаете, что это ошибка в одном, а не в другом, тогда я понимаю. - person user541686; 31.01.2014
comment
@Mehrdad: «Intel 64» - правильное название 64-битной архитектуры Intel; именно так он используется в документации Intel. Я предположил, что он может быть использован, поскольку вы упомянули разницу с использованием переключателя /Arch:IA32 или нет. Однако из этой документации я вижу, что переключатель назван неверно; он не выбирает архитектуру IA-32, а скорее выделяет отсутствие таких расширений, как SSE, SSE2 и AVX. В любом случае код компилятора для генерации инструкций для различных целевых функций должен отличаться. - person Eric Postpischil; 31.01.2014

arrow_upward
0
arrow_downward

Вы разыграли y1 как двойную, прежде чем снова разыграть ее на длинную. значение x не является значением «пола», а является округленным значением для этажа.

Та же логика применима к приведениям целых чисел и чисел с плавающей запятой. float x = (float) ((int) 1.5) даст другое значение для float x = 1.5

person James Malone    schedule 31.01.2014