Снова точность с плавающей запятой

Вчера я задал вопрос о том, почему я теряю точность арифметики с плавающей запятой. Я получил ответ о том, что это произошло из-за того, что промежуточные результаты хранятся в регистрах x87. Это было полезно, но некоторые детали все еще ускользают от меня. Вот вариант программы, которую я представил в предыдущем вопросе. Я использую VC ++ 2010 Express в режиме отладки.

int main()
{
    double x = 1.8939201459282359e-308; /* subnormal number */
    double tiny = 4.9406564584124654e-324; /* smallest IEEE double */
    double scale = 1.6;
    double temp = scale*tiny;
    printf("%23.16e\n", x + temp);
    printf("%23.16e\n", x + scale*tiny);
}

Это выводит

1.8939201459282369e-308
1.8939201459282364e-308

Первое значение соответствует стандарту IEEE. Присвоение переменной scale значения 2,0 дает правильное значение для обоих вычислений. Я понимаю, что temp в первом расчете является субнормальным значением и поэтому теряет точность. Я также понимаю, что значение scale*tiny хранится в регистре x87, который имеет больший диапазон экспоненты, и поэтому это значение имеет большую точность, чем temp. Чего я не понимаю, так это того, что при добавлении значения к x мы получаем правильный ответ из значения с более низкой точностью. Конечно, если более низкое значение точности может дать правильный ответ, то более высокое значение точности также должно дать правильный ответ? Это как-то связано с «двойным округлением»?

Заранее спасибо, это совершенно новая тема для меня, поэтому я немного борюсь.


c floating-point floating-accuracy ieee-754 floating-point-precision
person john    schedule 16.03.2013    source источник
comment
Следующее вполне может быть правдой, но для меня это совсем не очевидно: Конечно, если более низкое значение точности может дать правильный ответ, то более высокое значение точности также должно дать правильный ответ?   -  person NPE    schedule 16.03.2013
comment
На вашем месте я бы использовал long double в таких расчетах ...   -  person Rontogiannis Aristofanis    schedule 16.03.2013
comment
Как мы узнаем, что число с более низкой точностью не имеет случайного значения в последней цифре? Всегда есть 10% шанс попасть в ожидаемый.   -  person Bo Persson    schedule 16.03.2013
comment
@RondogiannisAristophanes Я хочу понять, что происходит.   -  person john    schedule 16.03.2013
comment
@BoPersson Ваш комментарий меня озадачивает, случайных цифр нет, все определяется. Плюс IEEE-754 с плавающей точкой является двоичной, а не десятичной.   -  person john    schedule 16.03.2013
comment
Я просто говорю, что если они оба ошибаются (= за пределами точности), один из них все равно может быть тем, которого вы ожидали.   -  person Bo Persson    schedule 16.03.2013
comment
Я ожидал результата IEEE-754. До недавнего времени у меня создалось впечатление, что моя платформа поддерживает IEEE-754 (но не по умолчанию).   -  person john    schedule 16.03.2013

Ответы (1)


arrow_upward
7
arrow_downward

Дело в том, что из-за большего диапазона экспонент эти два числа не являются субнормальными в представлении x87.

В представлении IEEE754

x    = 0.d9e66553db96f × 2^(-1022)
tiny = 0.0000000000001 × 2^(-1022)

но в представлении x87

x    = 1.b3cccaa7b72de × 2^(-1023)
tiny = 1.0000000000000 × 2^(-1074)

Теперь, когда 1.6*tiny вычисляется в представлении IEEE754, оно округляется до 0.0000000000002 × 2^(-1022), поскольку это наиболее близкое представимое число к математическому результату. Добавление этого к x приводит к

  0.d9e66553db96f × 2^(-1022)
+ 0.0000000000002 × 2^(-1022)
-----------------------------
  0.d9e66553db971 × 2^(-1022)

Но в представлении x87 1.6*tiny становится

1.999999999999a × 2^(-1074)

и когда это добавлено

  1.b3cccaa7b72de × 2^(-1023)
+ 0.0000000000003333333333334 × 2^(-1023)
-----------------------------------------
  1.b3cccaa7b72e1333333333334 × 2^(-1023)

результат, округленный до 53 значащих битов,

  1.b3cccaa7b72e1 × 2^(-1023)

с последним битом в мантиссе 1. Если это затем преобразовать в представление IEEE754 (где оно может иметь не более 52 бит в мантиссе, потому что это субнормальное число), поскольку оно находится ровно посередине между двумя соседними представляемыми числами 0.d9e66553db970 × 2^(-1022) и 0.d9e66553db971 × 2^(-1022) по умолчанию оно округляется до единицы с последним битом в нулевом значении.

Обратите внимание, что если бы FPU не был настроен для использования только 53 бита для значения, а полных 64 бита типа расширенной точности x87, результат сложения был бы ближе к результату 0.d9e66553db971 × 2^(-1022) IEEE754 и, следовательно, был бы округлен до этого.

Фактически, поскольку представление x87 имеет больший диапазон экспоненты, у вас есть больше битов для значений субнормальных чисел IEEE754, чем в представлении IEEE754, даже с ограниченным числом бит в мантиссе. Таким образом, результат вычисления здесь в x87 имеет на один значащий бит больше, чем в IEEE754.

person Daniel Fischer    schedule 16.03.2013
comment
Спасибо, Даниэль, рабочий пример был действительно тем, что мне было нужно. Итак, когда 1.b3cccaa7b72e1 × 2 ^ (- 1023) преобразуется обратно в IEEE-754, он округляется до 0.d9e66553db970 × 2 ^ (- 1022) вместо 0.d9e66553db971 × 2 ^ (- 1022)? Какой вообще режим округления для этой операции? - person john; 16.03.2013
comment
Верно. (Хотя я не знаю, округляется ли оно до IEEE754 для printf, printf также может использовать представление x87.) Режим округления по умолчанию в IEEE754 - это округление до четности, то есть последний бит значимого нуля . - person Daniel Fischer; 16.03.2013
comment
Привет, Даниэль, небольшое замечание: то, как вы описываете сложение в x87, рядом со словами «из-за ограничения значащих битов становится 0,0000000000003 × 2 ^ (- 1023)», звучит как сложение Cray (cs.nyu.edu/courses/fall03/G22.2420-001/lec4.pdf). Вместо этого x87 концептуально эквивалентен вычислению точной суммы (1.b3cccaa7b72e1333333333334 × 2 ^ (- 1023)) с последующим округлением. - person Pascal Cuoq; 16.03.2013
comment
@PascalCuoq Спасибо, не был уверен, как именно x87 работает в этой конфигурации. - person Daniel Fischer; 17.03.2013