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

Question

Вчера я спросил question о том, почему я терял точность в арифметике с плавающей запятой. Я получил ответ о том, как это произошло из-за того, что промежуточные результаты проводятся в регистре 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 мы получаем правильный ответ от значения меньшей точности. Разумеется, если значение более низкой точности может дать правильный ответ, то более точное значение точности также должно дать правильный ответ? Это как-то связано с «двойным округлением»?

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

Answer 1

Дело в том, что из-за большего диапазона экспоненциальности два числа не являются субнормальными в представлении 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, результат добавления будет ближе к результату IEEE754 0.d9e66553db971 × 2^(-1022) и, следовательно, округлен до этого ,

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