В двух словах
Вы говорите, что ваш компилятор — Visual C++ 2010 Express. У меня нет доступа к этому компилятору, но я понимаю, что он генерирует программы, изначально настраивающие ЦП x87 на использование 53-битной точности, чтобы как можно точнее эмулировать вычисления двойной точности IEEE 754.
К сожалению, «как можно ближе» не всегда достаточно близко. Исторические 80-битные регистры с плавающей запятой могут иметь значащую ограниченную ширину в целях эмуляции двойной точности, но они всегда сохраняют полный диапазон для экспоненты. Разница проявляется, в частности, при манипулировании денормалями (например, вашим y).
Что случается
Мое объяснение состоит в том, что в printf("%23.16e\n", 1.6*y); 1.6*y вычисляется как 80-битное число с уменьшенной значимостью и полной экспонентой (таким образом, это нормальное число), затем преобразуется в IEEE 754 с двойной точностью (что приводит к денормальному), а затем печатается.
С другой стороны, в printf("%23.16e\n", x + 1.6*y);, x + 1.6*y вычисляется со всеми 80-битными числами с пониженной значимостью и полной экспонентой (опять же, все промежуточные результаты являются нормальными числами), затем преобразуются в IEEE 754 с двойной точностью, а затем печатаются.
Это объясняет, почему 1.6*y печатается так же, как 2.0*y, но имеет другой эффект при добавлении к x. Напечатанное число является денормализованным с двойной точностью. Число, которое добавляется к x, представляет собой 80-битное нормальное число с уменьшенной значимостью и полной экспонентой (не то же самое).
Что происходит с другими компиляторами при генерации инструкций x87
Другие компиляторы, такие как GCC, не настраивают x87 FPU для работы с 53-битными мантиссандами. Это может иметь те же последствия (в этом случае x + 1.6*y будет вычисляться со всеми 80-битными числами с полной мантиссой и полной экспонентой, а затем преобразовываться в двойную точность для печати или сохранения в памяти). В этом случае проблема заметна даже чаще (не нужно задействовать денормали или бесконечные числа, чтобы заметить различия).
Эта статья Дэвида Моннио содержит всю необходимую информацию и многое другое.
Удаление нежелательного поведения
Чтобы избавиться от проблемы (если вы считаете ее таковой), найдите флаг, который указывает вашему компилятору генерировать инструкции SSE2 для операций с плавающей запятой. Они реализуют точно семантику IEEE 754 для одинарной и двойной точности.
person
Pascal Cuoq
schedule
15.03.2013
