Однако

[vak]

Знаете ли вы, что Интеловские процессоры врут при округлении?  Попробуйте умножить 3.00000000000000044 на 1.00000000000000022 - пентиум даёт 3.00000000000000089.  А между тем правильный ответ 3.00000000000000133.  Компилятор не виноват, он честно включает аппаратное округление к ближайшему числу.

То же самое при делении: 1 / 0.99999999999999989 должно давать 1.00000000000000022, а не 1.

Вывод: покупайте процессоры MIPS. :)

Upd: выяснилось, что баг проявляется только на Intel Xeon (W3520 и W3530) и Pentium 4.  Core Duo, Core2 Duo, i5 и AMD Opteron работают правильно.

Comments

[vak]

Да, проверил на Intel Core Duo - работает правильно.
А на Intel(R) Xeon(R) W3530 безбожно глючит.

[spamsink.livejournal.com]

Объясняю популярно, для невежд. Округление по умолчанию делается tie break to even. Если промежуточные вычисления делаются в long double, то 2*ε² оказываются за пределами разрядной сетки, 5*ε после нормализации выглядит как 2.5*ε, и округление делается вниз, к 2*ε.
Если промежуточные вычисления делаются в quad, то ничьей не оказывается, и получается честное округление к ближайшему, т.е. 3*ε.
И это не баг, а естественное ограничение конкретных реализаций FPU.

[vak]

В мипсе нет никаких quad, и тем не менее умножитель умудряется сообразить, что это не tie break. Забавно, что на разных интелах оно работает по разному.

[spamsink.livejournal.com]

Да, пардон, это не quad, это double-double (106 бит мантиссы). IEEE не специфицирует, с какой дополнительной точностью должно делаться умножение.

[lionet.livejournal.com]

Это не дабл-дабл, а extended precision (64 бит мантиссы).

[spamsink.livejournal.com]

Объясняю еще раз: 64 бит мантиссы недостаточно для понимания, что ε/2+2*ε² не является tie. Именно поэтому "старые" пентиумы этого и не понимают. А в новых анализируются все биты результата умножения двух double мантисс.

[lionet.livejournal.com]

double-мантисса — это 53 бита. Я же говорю про extended precision, которые 64-бита мантиссы. double-double нет в интелах.

[spamsink.livejournal.com]

Я вижу, имеется некоторое непонимание. Результат умножения в процессоре попадает в регистр, имеющий определенную длину мантиссы, которая больше, чем у double. Ненулевые разряды этого результата соответствуют 2*ε(double), ε(double)/2, и, если разрядность мантиссы у этого регистра вдвое больше, чем у double, &2*ε(double)². Если разрядность этого регистра такова, что &2*ε(double)² в нем нет (или в fpucw установлен соответствующий режим), то при сохранении регистра в память в формате double происходит округление по типу tie break to even, т.е. вниз, что мы и наблюдаем на старых версиях процессора. Если же регистр как минимум вдвое длиннее, чем мантисса double (и fpucw не запрещает анализировать его полностью), то происходит округление по типу round to nearest, т.е. вверх.

Говорил же, режимы нужно правильно ставить

[spamsink.livejournal.com]

void
set_fpu (unsigned int mode)
{
asm ("fldcw %0" : : "m" (*&mode));
}
int main(int ac, char **av) {
volatile double a = 3.00000000000000044;
volatile double b = 1.00000000000000022;
set_fpu(0x27F);
printf("%.17f\n", a * b);
}

Результат умножения будет немедленно, до попадания в регистр FPU округляться до double, и будет правильно на всех процессорах.

[vak]

Так-так, сейчас посмотрим, что за режимы...
Я пробовал fesetround(FE_TONEAREST) - не помогает.

[spamsink.livejournal.com]

По умолчанию 0x37F (3 - extended precision, 2 - double precision, 1 - reserved, 0 - single precision).

Я вот чего не пойму:

Xeon:
FPU mode 7f: 3.00000000000000000
FPU mode 17f: 3.00000000000000089
FPU mode 27f: 3.00000000000000133
FPU mode 37f: 3.00000000000000089

Core:
FPU mode 7f: 3.00000000000000133
FPU mode 17f: 3.00000000000000133
FPU mode 27f: 3.00000000000000133
FPU mode 37f: 3.00000000000000133

Что у них extended precision разная - понятно, но почему single precision так нагло игнорируется?

[vak]

Да, помогает. Почему же другие процессоры даже с 0x37F дают правильный результат?

[spamsink.livejournal.com]

Потому что их extended precision достаточна, чтобы eps*eps туда поместился.

[vak]

Я слышал раньше, что физики не любят Интел. Теперь понимаю, почему.

[vak]

На Оптероне эти биты тоже не влияют. Что же они, молча поменяли спецификацию такого важного регистра?

[spamsink.livejournal.com]

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