Неверное округление с плавающей запятой

В gcc 4.7.3 моя функция fegetround() возвращает FE_TONEAREST. Согласно справочнику по C++, это означает округление от нуля. По сути, это означает сохранение последнего бита, который был смещен при настройке точности мантиссы после умножения (поскольку она будет в два раза длиннее, чем должна быть). После этого сохраненный бит добавляется к окончательному результату мантиссы.

Например, умножение с плавающей запятой дает следующие результаты:

0x38b7aad5 * 0x38b7aad5 = 0x3203c5af

Мантисса после умножения

  1011 0111 1010 1010 1101 0101
x 1011 0111 1010 1010 1101 0101
-------------------------------
1[000 0011 1100 0101 1010 1110] [1]000 0101 1001 0101 0011 1001

Набор [23'b] содержит значащие цифры, тогда как набор [1'b] содержит последний сдвинутый бит. Обратите внимание, что мантисса для результата равна

[000 0011 1100 0101 1010 1111]

Последний бит переключился на 1, потому что набор [1'b1] был добавлен к объединенной мантиссе (набор [23'b]) из-за режима округления.

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

0x20922800 * 0x20922800 = 0x1a6e34c (check this on your machine)

  1010 0110 1110 0011 0100 1101
x 1010 0110 1110 0011 0100 1101
-------------------------------
01[01 0011 0111 0001 1010 0110 0][1]00 0000 0000 0000 0000 0000

Final Mantissas:       
Their Result:      01 0011 0111 0001 1010 0110 0
Correct Result(?): 01 0011 0111 0001 1010 0110 1

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


c++ floating-point floating-accuracy floating-point-precision
person Suedocode    schedule 06.12.2013    source источник
comment
«возвращает FE_TONEEAREST. Согласно справочнику по c++, это означает округление от нуля» -> вы уверены? Посмотрите еще раз на название.   -  person Pascal Cuoq    schedule 07.12.2013
comment
Я думаю, что всякий раз, когда возникает вопрос «Моя ошибка или ошибка Intel», вы можете дать ответ, не читая вопрос, и в 99% случаев вы будете правы.   -  person Kerrek SB    schedule 07.12.2013
comment
Если у вас нет Pentium первого поколения.   -  person dan04    schedule 07.12.2013
comment
Моя ошибка или интеловская? - Будь скромнее, мужик.   -  person    schedule 07.12.2013
comment
Хм, мой fegetround() возвращает 0, что, как я думал, где-то видел, означало FE_TONEAREST... Возможно, я ошибаюсь! РЕДАКТИРОВАТЬ: я только что успешно assert(fegetround() == FE_TONEAREST);ed.   -  person Suedocode    schedule 07.12.2013
comment
@ dan04: Сколько людей заполонили SO, чтобы спросить об их подразделениях с плавающей запятой на их новом Pentium, когда он был выпущен?   -  person Kerrek SB    schedule 07.12.2013
comment
@KerrekSB: существовало ли ТАК в 1994 году?   -  person Joe Z    schedule 07.12.2013
comment
@Aggieboy: округление до ближайшего означает округление от 0,5 до ближайшего четного числа. Кажется, это то, что происходит.   -  person Joe Z    schedule 07.12.2013
comment
@JoeZ: Его было не так просто найти через Lycos и Altavista, но он всегда занимал видное место в каталоге экспертных форумов Yahoo...   -  person Kerrek SB    schedule 07.12.2013
comment
@JoeZ В ссылке, которую я дал в OP, говорится, что округление x до ближайшего выбирает возможное значение, ближайшее к x, с промежуточными случаями округления от нуля. Я не вижу ничего упомянутого о четных числах.   -  person Suedocode    schedule 07.12.2013
comment
Ваш вопрос, вероятно, не был бы отклонен, если бы не заголовок.   -  person Reinstate Monica -- notmaynard    schedule 07.12.2013
comment
@Aggieboy: cppreference.com вводит в заблуждение. Режим округления IEEE-754 по умолчанию — округление до ближайшего, привязка к четному. См. здесь: en.wikipedia.org/wiki/   -  person Joe Z    schedule 07.12.2013
comment
@Aggieboy У вас есть основания верить этой отсылке?   -  person David Schwartz    schedule 07.12.2013
comment
@KerrekSB: Это было то же самое? На странице «о» StackExchange ( stackexchange.com/about ) говорится, что StackOverflow был создан в 2008 году. Или они просто ссылаются на SO как на отдельный сайт, а не форум на Yahoo?   -  person Joe Z    schedule 07.12.2013
comment
@JoeZ: Нет, я пошутил. SO действительно не так уж и стар. Но моя точка зрения заключалась в том, что, вероятно, большинство людей не сразу вышли и публично задались вопросом о делении с плавающей запятой, несмотря на ошибку. Это была досадная ошибка, но чтобы на нее наткнуться, потребовалась чрезвычайно специфическая и сложная работа, и, по-видимому, нашедший провел много тщательных исследований, чтобы убедиться, что это не его вина.   -  person Kerrek SB    schedule 07.12.2013
comment
@DavidSchwartz Я верил в это по той же причине, по которой Джо Зи верил в свою ссылку на вики. ^. ^   -  person Suedocode    schedule 07.12.2013
comment
@KerrekSB: Скорее, в 99,99999% случаев. Как минимум.   -  person Stephen Canon    schedule 07.12.2013

Ответы (2)


arrow_upward
6
arrow_downward

При округлении до ближайшего IEEE указывает, что округление округляется до четного. 0 четное, 1 нечетное, поэтому Intel прав.

person David Schwartz    schedule 06.12.2013
comment
Разве это не означает, что все неточные вычисления заканчиваются 0 в мантиссе? Почему не было первого? - person Suedocode; 07.12.2013
comment
@Aggieboy: округляется только от 1/2 до четного. Если вы даже на йоту выше 1/2, оно будет округлено, как в вашем первом примере. Ваш второй пример округлял точно 1/2. - person Joe Z; 07.12.2013
comment
Ага, это было объяснение, которое я искал! Я неправильно понял, что означает этот режим округления. - person Suedocode; 07.12.2013

arrow_upward
2
arrow_downward

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

Цитата по стандарту IEEE 754 (раздел 4.3.1):

roundTiesToEven, должно быть доставлено число с плавающей запятой, ближайшее к бесконечно точному результату; если два ближайших числа с плавающей запятой, заключающие в скобки непредставимый бесконечно точный результат, одинаково близки, то должно быть доставлено число с четной младшей значащей цифрой

В первом примере вы вычисляете квадрат (8.75794e-5), который (если представлен как 32-битное число с плавающей запятой) имеет следующий шестнадцатеричный шаблон: 0x38b7aad5.

Все 24 значащих бита (8.75794e-5):

0xb7aad5 = 1.0110111_10101010_11010101

Теперь после возведения в квадрат вы получите:

1.0000011_11000101_10101110_10000101_10010101_00111001

Примечательно, что в 99% случаев ваши вычисления будут выполняться на FPU (вероятно, x87), который работает с 80-битным форматом с плавающей запятой.

Руководство разработчика программного обеспечения для архитектур Intel® 64 и IA-32

(ПРОГРАММИРОВАНИЕ С ПОМОЩЬЮ X87 FPU):

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

Теперь вы хотите сохранить результат в 32-битном формате с плавающей запятой:

1.[0000011_11000101_10101110]10000101_10010101_00111001

и вот где режимы округления имеют значение. IEEE 754 определяет 4 из них, но давайте сосредоточимся на значении по умолчанию (с округлением до ближайшего (четного)), поскольку мы обсуждаем его здесь.

Теперь, когда ваш FPU имеет вычисленный результат (целый - у нас есть 64 значащих бита в 80-битном формате), он должен выполнить округление, чтобы число соответствовало 32 битам (24 значащих бита). Все 23 бита, которые должны быть явно сохранены, помещены в скобки выше.

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

1.[0000011_11000101_10101111]
and
1.[0000011_11000101_10101110]

но они ближе к

1.[0000011_11000101_10101111]

Вот почему значимость вашего результата равна 0x3203C5AF.

Теперь проблемный результат возведения в квадрат 2.4759832E-19 0x20922800.

24 значащих бита 2.4759832E-19:

0x922800 = 1.0010010_00101000_0000_0000

и в квадрате:

1.[0100110_11100011_01001100]10000000_00000000_0000000

И здесь действительно важна четная часть. Теперь ваше значение находится ровно посередине между:

1.[0100110_11100011_01001101]
and
1.[0100110_11100011_01001100]

Говорят, что более 2 значений заключают в скобки ваше значение. Из них теперь нужно выбрать хотя бы один (последний, так как lsb=0).

Теперь вы знаете, почему 24 бита вашего результата равны 0xA6E34C(ближайшие четные), а не 0xA6E34D(ближайшие, но нечетные)

person Artur    schedule 10.12.2013