Объясните 1 - 3 * (4/3 - 1) = 2,2204e-16 в Matlab

Вещественное число 4/3 не может быть представлено в формате двойной точности (или в любом другом двоичном формате с плавающей запятой), потому что оно не является двоично-рациональным. Таким образом, когда вы вычисляете 4/3 в MATLAB, полученное значение округляется до ближайшего представимого числа с двойной точностью, а именно:

1.3333333333333332593184650249895639717578887939453125

Вычитание 1 из этого значения является точным (известная теорема анализа ошибок FP гласит, что вычитание чисел, отличающихся друг от друга в два раза, является точным), поэтому результат 4/3 - 1 таков:

0.3333333333333332593184650249895639717578887939453125

Бывает, что результат умножения этого числа на три тоже точно представим:

0.9999999999999997779553950749686919152736663818359375

Наконец, вычитание из 1,0 также является точным (согласно теореме, на которую я ссылался ранее):

0.0000000000000002220446049250313080847263336181640625

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

Выполнение этого на научном калькуляторе - log₂(2.2204e-16) - дает (почти) ровно -52. Другими словами, Matlab хранит 52 бита точности в double, еще 5 бит для показателя степени (4 + знак) и один для знака мантиссы. Это соответствует реализации IEEE 754: 53 бита для мантиссы (52 + знак) и 5 ​​для экспоненты. Все хорошо! Как всегда, вы должны проверить, достаточно ли близки два числа с плавающей запятой, а не точно ли они равны. Соответствующий пример в Matlab будет выглядеть так:

if abs(x - y) < 1e-15
      % Some code to execute when x and y are approximately equal
else
      % Some other code
end

Заявление о соответствии стандарту IEEE 754 взято из статьи в Википедии.

Проблемы начинаются с вычисления 4/3. Точный ответ не может быть представлен ни конечным числом десятичных цифр, ни конечным числом битов. Результат будет сохранен как 4/3+r, где r — небольшое число со знаком абсолютного значения, представляющее разницу между реальным значением 4/3 и ближайшим 64-битным двоичным числом с плавающей запятой IEEE 754 к 4/3, ошибке округления.

Вычитание 1 дает 1/3+r. Умножение на 3 дает 1+3r. Вычитание из 1 дает -3r.

Конечный результат в -3 раза превышает ошибку округления в исходном представлении 4/3.