Вздумалось кому-то (не мне) проверять, что матрица между вычислениями не поменялась. Проверять решил просто: считал определитель, сохранял значение и в нужный для проверки момент вычислял определитель опять. Если определитель не изменился, то можно спать спокойно.
На этом математика кончается и начинается песня. Код попадает ко мне — и начинаются глюки на самой свободной ОС, допиленной сумрачным нордическим гением аж до зелёного хамелеона, одиннадцатой версии и второго сервис-пака.
В результате отладки дохожу до такого кода:
double a = det(M);
assert(a == det(M));
Ассерт срабатывает. Ладно, добавляю строчку:
assert(det(M) == det(M));
Ассерт не срабатывает. Функция всегда возвращает одно и то же значение. Добавляю:
double diff = a - det(M);
Результат равен нулю. Причём строго нулю, посмотрел побайтово. Та-ак… Похоже, что имеем вещественное число, в общем случае не равное самому себе. Уже интересно…
double a = det(M);
double b = det(M);
assert(a == b);
Ассерт не срабатывает. Пора в дурку…
Ларчик открывался просто. В сопроцессоре все числа обрабатываются в 10-байтовом формате, а double, как известно, 8 байт. Разработчики самого безглючного компилятора возвращали значение в голове стека сопроцессора и забыли нормализовать его до 8 байт. Нормализация происходила только в случае сохранения значения в переменной. Хвост в 2 байта добавлял несколько знаков к мантиссе и вызывал все эти спецэффекты.