Множество вещественных чисел всюду плотное. Это значит, что на любом конечном интервале таких чисел будет бесконечное количество (и эта бесконечность будет даже не счётной, а мощности континуума). Поэтому, в отличие от целых чисел, невозможно придумать такой способ представления вещественного числа в памяти компьютера, который позволял бы представить все вещественные числа даже на небольшом интервале. Любое представление даст лишь конечное количество чисел на интервале, а все остальные придётся округлять до одного из этих значений, т.е. хранить с погрешностью. Хороший способ представления должен обеспечивать не слишком большую погрешность, но при этом позволять хранить числа в большом диапазоне, требовать не слишком много памяти, а ещё желательно, чтобы относительная погрешность чисел в разных частях диапазона была примерно одинаковой.
Это не такая простая задача, и первые попытки её решить были не очень успешными. Работа с вещественными числами в некоторых старых системах требовала глубокого знания способа их хранения, чтобы случайно не выполнить какую-нибудь внешне безобидную операцию, которая приведёт к некорректному результату или к ошибке. Но опыт был постепенно накоплен, и в 1985 году был принят стандарт IEEE 754, описывающий способ хранения вещественных чисел и операции с ними. Этот стандарт быстро завоевал популярность и с некоторыми дополнениями используется в большинстве современных систем.
Стандарт IEEE 754 настолько удачен, что позволяет в большинстве случаев работать с вещественными числами, особо не задумываясь об их внутреннем устройстве, результат всё равно будет таким, какой мы ожидаем. Тем не менее, существуют ситуации, когда вычисления, выполняемые в рамках стандарта, дают результаты, далёкие от интуитивно ожидаемых. А так как современный программист, в отличие от своих предшественников эпохи до IEEE 754, обычно плохо понимает, как выполняются вещественные операции на двоичном уровне, такие ситуации его полностью обескураживают и часто даже списываются на какие-то ошибки в работе процессора, системы и т.п. Хотя никакие это не ошибки, а прямое следствие особенностей стандарта IEEE 754.
В этой статье я собрал коллекцию таких ситуаций, которые возникают на платформах .NET и .NET Framework в Windows на процессорах семейства x86 (включая 64-разрядные версии). Обе этих платформы используют IEEE 754 для представления вещественных чисел. Большинство этих ситуаций воспроизведутся и на других программных и аппаратных платформах, использующих IEEE 754, но своя специфика у выбранных мною платформ есть. По каждой ситуации будут объяснения, почему результат получается именно такой.
Сразу сообщаю неприятную новость: вещественные вычисления в .NET и .NET Framework реализованы разными способами, и в некоторых случаях мы будем сталкиваться с тем, что один и тот же код на этих платформах даёт разные результаты. О причинах этого мы тоже будем говорить. Более того, в .NET Framework по-разному реализуются вычисления при работе на 32-разрядных и 64-разрядных процессорах, так что мало того, что ваша программа может начать работать по-другому при переходе с .NET Framework на .NET, но даже при переходе от 32-разрядной системы на 64-разрядную вы можете столкнуться с тем же самым. И, что самое неприятное, 32-разрядные приложения для .NET Framework иногда дают разные результаты в конфигурациях Debug и Release.
В примерах к этой статье под .NET Framework я буду подразумевать то, что код откомпилирован для платформы x86. К вопросу о том, что происходит в случае x64 и Any CPU, я вернусь в конце статьи. К статье прилагаются два проекта: FPNumbers содержит код всех примеров для .NET, а FPNumbersFW – всех примеров для .NET Framework. Я использовал версии .NET 6 и .NET Framework 4.8, но, насколько я знаю, в других версиях (за исключением, возможно, самых ранних версий .NET Framework) результаты будут такие же.
Сравнение значений различных типов
В своём базовом варианте стандарт IEEE 754 содержит описание двух форматов хранения вещественных чисел: число одинарной точности (размер 4 байта, минимальное отличное от нуля значение 1.4e-45, максимальное значение 3.4e38, точность 6-9 десятичных знаков) и число двойной точности (размер 8 байт, минимальное отличное от нуля число 4.9e-324, максимальное число 1.8e308, точность 15-17 десятичных знаков). В .NET (Framework) им соответствуют типы System.Single (в C# для него имеется синоним float) и System.Double (синоним double).
Начнём с простого вопроса: что будет, если сравнить одинаковые вещественные числа, хранящиеся в переменных различного типа?
Console.WriteLine("Equality of double and float");
f = 0.109344482421875f;
d = 0.109344482421875;
Console.WriteLine($"f == d for 0.109344482421875: {f == d}");
f = 0.1f;
d = 0.1;
Console.WriteLine($"f == d for 0.1: {f == d}");
Результат в .NET:
Equality of double and float
f == d for 0.109344482421875: True
f == d for 0.1: False
Результат в .NET Framework:
Equality of double and float
f == d for 0.109344482421875: True
f == d for 0.1: False
Здесь и в остальных примерах мы будем придерживаться соглашения, что переменные с названиями на f имеют тип float, на d – double.
Обе платформы показывают одинаковый результат, и этот результат прямо противоположен тому, который ожидается интуитивно. Давайте разберём подробно, какие результаты мог бы ожидать человек, знакомый с вещественными числами только поверхностно. Точность значения типа float составляет 6-9 десятичных знаков, а типа double – 15-17. В числе 0.109344482421875 мы имеем пятнадцать знаков, поэтому при записи в переменную типа double это значение сохранится без изменений, а при записи в переменную типа float будут отброшены как минимум шесть знаков, и в лучшем случае будет записано значение 0.109344482. При сравнении же значение типа float будет расширено до double, недостающие разряды будут заполнены нулями, и будут сравниваться числа 0.10934448200000 и 0.109344482421875, в результате мы получим false. А число 0.1 будет записано в обе переменные без усечений, поэтому и при сравнении получим, что значения равны.
Чтобы понять, почему такие простые и понятные рассуждения дают неверный результат, надо погрузиться в теорию и разобраться с тем, как устроены вещественные числа на бинарном уровне.
Для начала вспомним, что такое вообще число с плавающей десятичной точкой. Это специальная форма записи вещественного числа, состоящая из мантиссы и экспоненты. Мантисса – это число, модуль которого лежит в интервале 1 (включительно) до 10 (не включительно). Экспонента – это целое число, показывающее, на десять в какой степени нужно умножить мантиссу, чтобы получить требуемое число. В текстовом представлении такого числа записывается сначала мантисса, потом экспонента, они разделяются латинской буквой E (допускается любой регистр). Например, 10 – это 1e1, 555 – 5.55e2, 0.00034 – 3.4e-4 и т.п. Точка называется плавающей, потому что в мантиссе она всегда стоит после первого знака, но как бы «уплывает» вправо или влево на количество знаков, указанное в экспоненте.
Компиляторы обычно не ограничивают диапазон, в котором должна лежать мантисса. Так, например, числа 55.5e1 и 0.034e-2 будут правильно распознаны как 5.55e2 и 3.4e-4, но такая запись считается не совсем корректной. Корректная запись позволяет легко понять, сколько десятичных знаков содержит число, а это принципиально там, где числа записываются с какой-то определённой точностью.
Далеко не все вещественные числа могут быть записаны в такой форме без потери точности. Во-первых, любое число, записанное в форме с плавающей десятичной точкой, будет рациональным, потому что в его основе лежит рациональная десятичная дробь. Все иррациональные числа могут быть переданы только приближённо. Во-вторых, из рациональных чисел (т.е. чисел, представимых в виде результата деления целого числа на другое целое число) в виде десятичной дроби могут быть переданы только те, знаменатель которого равен 2^N*5^M, где N и M – некоторые целые неотрицательные числа. Числа же вида 1/3, 5/7 и т.п. не могут быть записаны конечной десятичной дробью. В этом случае получается периодическая дробь, т.е. такая десятичная дробь, в которой последовательность из одной или нескольких цифр повторяется бесконечное число раз. Кроме того, на практике количество разрядов в записи вещественного числа с плавающей десятичной точкой обычно ограничивается из каких-либо соображений. Так, например, если у нас три десятичных разряда, то число 12345 будет записано как 1.23e4, т.е. тоже приближённо.
Вещественные числа в формате стандарта IEEE 754 – это тоже числа с плавающей точкой, только не десятичной, а двоичной. Их мантисса и экспонента записываются в двоичном виде. Из 32 бит числа одинарной точности самый старший бит содержит знак числа (0 – положительное, 1 – отрицательное), затем идут 8 бит экспоненты, а младшие 23 бита содержат мантиссу. Число двойной точности устроено так же, только в нём на экспоненту отводится 11 бит, а на мантиссу 52 бита.
Двоичная дробь, хранящаяся в мантиссе, похожа на десятичную, только каждый её разряд после точки отвечает за соответствующую степень двойки: 1/2, 1/4, 1/8 и т.д. Т.е., например, 0.1 – это в десятичном виде 0.5, 0.01 – это 0.25, 0.0101 – 0.3125 (1/4 + 1/16) и т.д.
Как и в случае с десятичной плавающей точкой, мантисса должна лежать в диапазоне от 1 до 10, только теперь это двоичные 1 и 10 (т.е. 1 и 2 в десятичном представлении). Это значит, что у мантиссы тот разряд, который идёт до точки, всегда равен 1. А раз так, то можно не тратить бит на хранение этой единицы, а хранить в мантиссе только те разряды, которые стоят после точки (это называется «нормализованная запись». Это даёт нам один дополнительный бит в дробной части мантиссы, т.е. повышает точность представления числа в два раза.
Точность мантиссы фиксирована, но только в смысле количества двоичных знаков. Упоминаемые значения 6-9 знаков для чисел одинарной точности и 15-17 для двойной – это, как мы увидим дальше, достаточно условные диапазоны, в некоторых случаях точность в десятичном выражении может быть как выше, так и ниже этих значений.
Экспонента – это целое число, которое может быть как положительным, так и отрицательным, но представление экспоненты отличается от обычного представления знаковых целых чисел. Здесь используется т.н. целое со смещением, т.е. чтобы получить значение экспоненты, надо сначала интерпретировать эту комбинацию бит как обычное беззнаковое целое, а потом вычесть из него смещение. Для чисел одинарной точности смещение равно 127, для чисел двойной точности – 1023. Это даёт нам диапазон от -127 до 128 для чисел одинарной точности и от -1023 до 1024 для чисел двойной точности. Но максимальные значения экспоненты зарезервированы для особого применения (об этом чуть позже), поэтому на диапазоны получаются от -127 до 127 и от -1023 до 1023.
Давайте грубо оценим, какие минимальные и максимальные значения могут храниться в таких форматах. Максимальное значение мантиссы получается, когда она содержит все единицы. С учётом мысленно добавленной единицы перед точкой получаем значение, которое можно оценить сверху как 2. Соответственно, максимальное значение для числа одинарной точности примерно равно 2*2^127=3.4e38, двойной точности – 2*2^1023=1.8e308. Минимальное значение мантиссы равно 1 (оно достигается, когда все биты равны нулю – это означает равенство нулю дробной части, к которой мысленно добавляется единица перед точкой). Соответственно, получаем 2^-127=5.9e-39 для одинарной точности и 2^-1023=1.1e-308 для двойной.
Вернувшись к параметрам типов, приведённым на начале раздела, мы увидим, что верхние пределы мы вычислили правильно, а нижние у нас оказались значительно завышены. Значит, что-то в своих рассуждениях мы не учли. А не учли мы то, что на самом деле при минимально возможном значении экспоненты для мантиссы используется денормализованная запись, т.е. мысленно добавляемая к ней целая часть равна не 1, а 0. Это позволяет расширить диапазон за счёт снижения точности представления самых маленьких значений.
Денормализованная запись позволяет естественным образом хранить нулевое значение. Когда все биты мантиссы и экспоненты равны нулю, для числа одинарной точности получается значение 0*2^-127, т.е. 0. Никакие другие комбинации значений мантиссы и экспоненты не могут дать 0, потому что при любом другом значении экспоненты все нули в мантиссе будут означать 1, а не 0 за счёт дополнительного подразумеваемого целого разряда.
Таким образом, минимально возможное по модулю отличное от нуля число получается тогда, когда экспонента имеет своё минимальное значение, а в мантиссе только самый младший бит содержит 1. Для чисел одинарной точности младший бит в мантиссе соответствует 2^-22, а минимальное значение, соответственно, равно 2^-(22+127) = 1.4e-45. Аналогично, минимальное ненулевое число для значения двойной точности равно 2^-(51+1023) = 4.9e-324. Теперь наши теоретические оценки соответствуют тому реальным параметрам типов.
Вернёмся к максимальным значениям экспоненты, зарезервированным для особых случаев. Если все биты мантиссы равны нулю, а экспонента имеет своё максимальное значение, то такое значение называется бесконечностью. Если при вычислении получается значение, которое превышает максимально возможное для данного типа, оно заменяется на бесконечность. Любое другое значение мантиссы при максимальной экспоненте называется NaN – Not a Number, не число. Это значение является результатом операций, которые невозможно выполнить, например, взятия квадратного корня или логарифма от отрицательного числа. Константы float.NaN и double.NaN содержат мантиссу, в которой только старший бит равен 1, остальные нули, а также единице равен бит знака (для значения NaN он игнорируется). Но это не единственно возможная комбинация бит, которая интерпретируется как NaN. Если вы используете какие-то сторонние библиотеки для вычислений, особенно нативные, вы можете получит другой NaN, так что сравнивать эти значения побитно бессмысленно. Более того, по стандарту IEEE 754 любое сравнение, в котором участвует NaN, должно давать false, даже если это сравнение значения с самим собой. Проиллюстрируем это требование простым примером.
Console.WriteLine("NaN comparison");
f = float.NaN;
Console.WriteLine($"f == f: {f == f}");
Console.WriteLine($"f != f: {f != f}");
d = double.NaN;
Console.WriteLine($"d == d: { d == d}");
Console.WriteLine($"d != d: {d != d}");
Результат в .NET:
NaN comparison
f == f: False
f != f: True
d == d: False
d != d: True
Результат в .NET Framework:
NaN comparison
f == f: False
f != f: True
d == d: False
d != d: True
Есть даже рекомендация проверки того, является ли значение переменной NaN – сравнить переменную с самой собой. Если результат сравнения отрицательный, то это NaN. Впрочем, в .NET и .NET Framework есть методы float.IsNaN и double.IsNaN, которые позволяют выполнить эту проверку более наглядным способом.
Согласно первоначальной версии стандарта IEEE 754 любые вычисления с NaN в
качестве одного из аргументов должны давать в результате NaN. Но многие реализации вещественной арифметики при возведении 1 в степень не делали проверку на NaN и в нарушение требований стандарта возвращали 1. Поэтому в версии стандарта от 2008 для этой операции было сделано исключение, она теперь должна возвращать единицу. Так как .NET Framework появился до 2008 года, а .NET – после, они ведут себя в этой ситуации по-разному.
Console.WriteLine("1 to power of NaN");
Console.WriteLine($"1^NaN={Math.Pow(1, double.NaN)}");
Результат в .NET:
1 to power of NaN
1^NaN=1
Результат в .NET Framework:
1 to power of NaN
1^NaN=NaN
И последний момент, который осталось разобрать, – это знак числа. Для обычных чисел это очевидная вещь, но формат позволяет иметь знак даже таким значениям, как 0 и бесконечность. Бесконечность, как я уже говорил, получается в тех случаях, когда результат операции превышает максимально возможное значение. Такое может получиться, например, при перемножении двух больших чисел. Если эти числа имеют одинаковый знак, то получится плюс бесконечность, если разный – минус бесконечность. Аналогичным образом ноль может получиться, если результат операции, в принципе, не должен быть нулём, но его модуль меньше, чем минимальное отличное от нуля значение. В этом случае может получиться минус ноль, если понятно, что результат должен быть отрицательным. Это иллюстрируется следующим примером:
Console.WriteLine("Plus and minus");
d = 1e200;
d *= 1e200;
Console.WriteLine($"1e200 * 1e200 = {d}");
d = -1e200;
d *= 1e200;
Console.WriteLine($"-1e200 * 1e200 = {d}");
d = 1e-200;
d *= 1e-200;
Console.WriteLine($"1e-200 * 1e-200 = {d}");
d = -1e-200;
d *= 1e-200;
Console.WriteLine($"-1e-200 * 1e-200 = {d}");
Console.WriteLine($"-0 == 0: {d == 0}");
Результат в .NET:
Plus and minus
1e200 * 1e200 = ∞
-1e200 * 1e200 = -∞
1e-200 * 1e-200 = 0
-1e-200 * 1e-200 = -0
-0 == 0: True
Результат в .NET Framework:
Plus and minus
1e200 * 1e200 = ∞
-1e200 * 1e200 = -∞
1e-200 * 1e-200 = 0
-1e-200 * 1e-200 = 0
-0 == 0: True
Примечание: в шрифте, который используется консолью в Windows по умолчанию, отсутствует символ «∞», поэтому вместо него вы можете увидеть вопросительный знак.
Тут может показаться, что .NET Framework, в отличие от .NET, даёт в последней операции положительный ноль вместо отрицательного. На самом деле там тоже отрицательный ноль, а различия заключаются только в реализации метода double.ToString – видимо, разработчики .NET Framework просто забыли учесть этот момент.
Последняя операция сравнения добавлена в пример для того, чтобы показать, что, хотя -0 и 0 различимы, их сравнение даёт интуитивно ожидаемый положительный результат.
Теперь мы можем вернуться к нашему первому примеру и разобраться, почему он даёт именно такой результат. Число 0.109344482421875 – это 3583/32768, а 32768 – это 2^15, т.е. для хранения такой дроби требуется всего пятнадцать двоичных знаков после точки. Тип float имеет таких знаков 23, поэтому указанное число даже в переменную типа float записывается точно, без округления, несмотря на то, что в десятичной записи оно имеет 15 знаков. При расширении этого числа до двойной точности младшие разряды мантиссы заполняются нулями. Но в представлении этого числа с двойной точностью там и должны стоять нули, поэтому сравнение даёт положительный результат.
Число 0.1 – это 1/(2*5), а 5 не является степенью двойки. Соответственно, при попытке записать это число в виде двоичной дроби получаем бесконечную периодическую дробь. Она будет равна 0.0(0011). Соответственно, при записи такой дроби в переменную типа float бесконечная дробь обрежется на 23 знаках, в переменную типа double – на 52 знаках, и избыточные (по сравнению с float) биты мантиссы в значении типа double не все будут содержать нули. А вот при расширении значения типа float до double эти разряды будут заполнены нулями, и мы получим не совсем такое число с закономерным отрицательным результатом сравнения.
Для интереса посмотрим, насколько значение 0.1f отличается от 0.1d.
Console.WriteLine("Float and double 0.1");
d = 0.1f;
Console.WriteLine($"0.1f = {d}");
d2 = 0.1;
Console.WriteLine($"0.1d = {d2}");
Console.WriteLine($"0.1f > 0.1d: {d > d2}");
Результат в .NET:
Float and double 0.1
0.1f = 0.10000000149011612
0.1d = 0.1
0.1f > 0.1d: True
Результат в .NET Framework:
Float and double 0.1
0.1f = 0.100000001490116
0.1d = 0.1
0.1f > 0.1d: True
На первый взгляд кажется, что 0.1f в .NET и .NET Framework переводится в значение двойной точности немного по-разному. На самом деле это не так: значения там абсолютно одинаковые, а по-разному реализован метод double.ToString: в .NET он обычно выдаёт на один-два знака больше, чем в .NET Framework. Мы будем регулярно сталкиваться с этим и в дальнейших примерах, поэтому я больше не буду заострять на этом внимание.
Значение 0.1d выводится как 0.1, хотя мы уже знаем, что там не 0.1, а та часть этого числа, которая может уместиться в пятидесяти двух разрядах. Но метод ToString достаточно интеллектуален, чтобы вывести это значение в удобном для пользователя виде.
Обратите внимание, что 0.1f больше, чем 0.1d. На первый взгляд это кажется невозможным: когда число одинарной точности расширяется до двойной точности, то дополнительные разряды мантиссы просто заполняются нулями, а это не может привести к увеличению числа, поэтому 0.1d, у которого эти разряды не все равны нулю, казалось бы, должно быть больше. Само по себе это рассуждение правильное, но оно не учитывает тот факт, что, когда число 0.1 превращается в число одинарной точности, там не просто отсекаются лишние разряды, а выполняется округление к ближайшему числу, и в данном случае это округление вверх. Если мы посмотрим два ближайших к 0.1 числа, которые могут быть без потерь записаны в виде числа одинарной точности, то это будет примерно 0.10000000149 и 0.09999999404. Для первого из них погрешность будет примерно 1.5e-9, для второго – 6e-9. Поэтому выбирается то число, которое больше, так как это минимизирует погрешность. Таким образом, увеличение происходит не тогда, когда 0.1f расширяется до двойной точности, а когда 0.1 превращается в 0.1f.
Связь с целыми числами и порядок байт
Вещественные числа в формате IEEE 754 имеют одну интересную особенность: если соответствующую комбинацию байт увеличить или уменьшить на единицу по правилам работы с целыми числами, то мы получим следующее или предыдущее вещественное число. Разберём это более подробно.
Начнём с комбинации, когда все биты числа равны нулю. В формате IEEE 754 это означает положительный ноль. Сделаем инкремент этой комбинации, рассматривая её как целочисленное значение, получим комбинацию, в которой самый младший бит равен 1, остальные нули. Интерпретируя это как вещественное число, получим значение с минимально возможной экспонентой и мантиссой с единицей в младшем разряде. Это будет минимально возможное отличное от нуля число в данной записи. Сделав ещё один целочисленный инкремент, получим число с мантиссой 0…010, т.е. следующее возможное число, затем – с мантиссой 0…011 и т.д. В конце концов дойдём до ситуации, когда все биты мантиссы окажутся равны единицы. Следующий инкремент приведёт к тому, что все биты мантиссы обнулятся, а единица окажется в младшем бите экспоненты. В этом случае нулевая мантисса трактуется как 1, и мы получаем минимально возможное нормализованное число. Продолжая операции инкремента, мы снова дойдём до ситуации, когда все биты мантиссы равны 1. Следующий инкремент приведёт к тому, что мантисса снова обнулится (другими словами, станет равной 1), а экспонента увеличится на 1, что опять даст нам следующее возможное число. Так будет продолжаться, пока мы не дойдём до значения экспоненты, которое на 1 меньше максимально возможного, а все биты мантиссы равны 1. Это будет максимально возможное в данной записи число. Следующий инкремент обнулит мантиссу, доведёт значение экспоненты до максимально возможного, и мы получим положительную бесконечность. Дальнейшие инкременты будут увеличивать мантиссу, и мы будем получать различные варианты «положительных» NaN (я взял слово «положительный» в кавычки потому, что для NaN знаковый бит не имеет смысла, но я хотел подчеркнуть, что на данном этапе мы будем получать только те варианты NaN, у которых этот бит равен 0). Наконец, мы дойдём до комбинации, при которых и мантисса, и экспонента состоят из одних единиц – это будет последний «положительный» NaN. После следующего инкремента единица будет только в самом старшем бите, отвечающем за знак, а все биты мантиссы и экспоненты обнулятся. Это будет отрицательный ноль. На следующем шаге мы получим число с единицами в знаковом разряде и в младшем разряде мантиссы, т.е. наименьшее по модулю отрицательное число, отличающееся от нуля. Продолжая инкременты, будем получать всё большие по модулю отрицательные числа, пока не дойдём до отрицательной бесконечности, после которой пойдут «отрицательные» NaN. Последний «отрицательный» NaN, в котором все биты равны 1, после инкремента даст положительный ноль, с которого мы и начали.
Описанная выше закономерность – это не забавное случайное совпадение, а запланированная и документированная особенность стандарта IEEE 754, которой вы можете пользоваться в своих проектах. Кроме всего прочего, это означает, что порядок байт в вещественном числе привязан к порядку байт в целом числе на данной платформе. Память адресуется побайтно, если число занимает несколько байт, то оно хранится в нескольких последовательных ячейках, каждая из которых имеет свой адрес. Если байты хранятся в порядке от младшего к старшему, это называется Big Endian (BE), если от старшего к младшему, то это Little Endian (LE) (я не нашёл какого-то общепринятого перевода этих терминов на русский язык, поэтому дальше буду пользоваться английскими аббревиатурами). Рассмотрим это на примере двухбайтного числа 12345 (0x3039 в шестнадцатеричном представлении). При использовании BE байты будут иди в таком порядке: 39 30. При использовании LE – в противоположном порядке: 30 39.
Процессоры архитектуры x86 (включая x86-64) используют LE. Это значит, что, например, четырёхбайтное вещественное число одинарной точности будет храниться следующим образом. Байт с нулевым смещением хранит восемь младших бит мантиссы. Байт со смещением 1 хранит следующие 8 бит мантиссы. В младших семи битах байта со смещением 2 хранятся старшие семь бит мантиссы, а в его старшем бите – младший бит экспоненты. В байте со смещением 3 семь младших бит хранят семь старших бит экспоненты, а старший бит хранит знак. Восьмибайтное число двойной точности хранится аналогичным образом.
Когда закладывались принципы обмена данными через сеть, процессоры с LE были менее распространены, чем сейчас, поэтому при передаче данных через сеть был принят вариант BE. Таким образом, если программа, выполняемая на процессоре x86, получает по сети бинарные данные, содержащие вещественные числа, и протокол передачи этих данных следует сетевым правилам, то байты в принятых числах будут идти в обратном порядке, и их надо будет переставить, прежде чем будет получено корректное представление числа. Впрочем, в .NET (Framework) редко приходится делать это самостоятельно: порядок байт – это вопрос протокола уровня представлений (6-й уровень в модели OSI), а код .NET обычно использует более высокоуровневые библиотеки, которые сами заботятся о перестановке байт, если это требуется. Но если вы отправляете или получаете данные непосредственно через сокеты, этот момент нужно учитывать. Другой случай, когда вам может потребоваться переставить байты в вещественном числе – если приходится разбирать бинарный файл, сформированный в системе, которая использует BE.
Умножение 0.1 на 10
Итак, при использовании типов float и double число 0.1 может быть представлено только приближённо, с некоторой погрешностью. А что будет, если мы это число умножим на 10? Поучится ли в результате 1? Чтобы узнать возможную погрешность как можно точнее, я вычту результат этого умножения из единицы.
Console.WriteLine("0.1*10");
f = 0.1f;
f = 1 - 10 * f;
Console.WriteLine($"1 - 0.1f*10: {f}");
d = 0.1;
d = 1 - 10 * d;
Console.WriteLine($"1 - 0.1d*10: {d}");
Результат в .NET:
0.1*10
1 - 0.1f*10: 0
1 - 0.1d*10: 0
Результат в .NET Framework:
0.1*10
1 - 0.1f*10: -1.490116E-08
1 - 0.1d*10: 0
В операциях с числами двойной точности мы в обоих случаях получили правильное значение 0, а при использовании чисел одинарной точности правильный результат показал только .NET.
Чтобы понять, откуда взялась разница между .NET и .NET Framework, поставим точку останова и посмотрим, какой код выполняет процессор в каждом из этих случаев.
Код в .NET:
f = 0.1f;
vmovss xmm0,dword ptr [Base+0ED0h]
vmovss dword ptr [rbp+12Ch],xmm0
f = 1 - 10 * f;
vmovss xmm0,dword ptr [rbp+12Ch]
vmulss xmm0,xmm0,dword ptr [Base+0F18h]
vmovss xmm1,dword ptr [Base+0EB0h]
vsubss xmm0,xmm1,xmm0
vmovss dword ptr [rbp+12Ch],xmm0
Код в .NET Framework:
f = 0.1f;
mov dword ptr [ebp-40h],3DCCCCCDh
f = 1 - 10 * f;
fld dword ptr [ebp-40h]
fmul dword ptr [Base+0A00h]
fld1
fsubrp st(1),st
fstp dword ptr [ebp-40h]
Здесь я для краткости заменил словом Base базовый адрес, от которого отсчитывается место хранения констант в коде. В реальности оно отсчитывается от места расположения кода метода FPNumbers.Program.Main.
Даже если вы мало что понимаете в ассемблере, вам должно быть заметно, что код .NET совсем не похож на код .NET Framework. Чтобы понять, почему так, надо разобраться, какие операции с числами IEEE 754 поддерживают процессоры x86.
Исторически процессоры Intel не умели самостоятельно работать с вещественными числами, но могли отдавать эту работу «на аутсорс» другой микросхеме, которая называлась математическим сопроцессором, или просто сопроцессором. Сопроцессор был необязательным элементом системы, в случае его отсутствия попытка использования команд для работы с вещественными числами вызывала прерывание (в данном случае прерывание можно рассматривать как некий низкоуровневый аналог исключения).
В первой половине 1990-х годов Intel добавила в свои процессоры FPU (Floating Point Unit), который выполнял все те же команды, что и ранее сопроцессор. С тех пор все процессоры Intel и совместимые с ними поддерживают операции с вещественными числами. FPU до сих пор иногда называют сопроцессором, хотя это не совсем корректно.
FPU имеет стековую систему команд, стек содержит 10-байтные вещественные числа. Этот формат не является частью стандарта IEEE 754, но устроен по принципам, близким к описанным в этом стандарте. На мантиссу в нём отводится 64 бита, на экспоненту 15 бит. Но, в отличие чисел одинарной и двойной точности, здесь явно хранится единица перед точкой (или ноль для денормализованных чисел), поэтому фактически такие числа имеют точность 63 двоичных знака после точки. Это даёт диапазон от 3.6e-4951 до 1.18e4932 и ориентировочную точность 19-20 десятичных знаков.
Все вычисления FPU выполняет над этими десятибайтными регистрами. Совместимость со стандартом IEEE 754 заключается в том, что у FPU есть команды, которые загружают в его регистры числа, представленные в формате одинарной или двойной точности из памяти и выгружают числа из его регистров, преобразуя в значения одинарной или двойной точности. При загрузке значений в регистры избыточные биты мантиссы дополняются нулями, при выгрузке мантисса округляется до точности целевого типа, а если число в регистре слишком велико или слишком мало для целевого типа, оно заменяется на бесконечность или ноль.
Существует возможность также загружать и выгружать десятибайтные значения без потерь, т.е. реализовывать полноценную работу с типом, более точным, чем double. Некоторые компиляторы имеют специальный тип данных для этого, например, long double в компиляторах C/C++ от Intel или Extended в Delphi. Но компиляторы, ориентированные на мультиплатформенность, обычно избегают этого типа, так как он не является частью какого-либо стандарта и поддерживается не всеми аппаратными платформами (у меня нет информации о том, чтобы такой формат поддерживался где-либо, кроме x86). Ни .NET Framework, ни .NET не имеют типов для работы с 10-байтным типом FPU.
В 1999 году процессоры Intel получили дополнительный набор инструкций, известный как SSE. Эти инструкции ориентированы на выполнение вычислений в задачах, связанных с мультимедиа. Для таких задач характерна ситуация, когда одну и ту же операцию надо повторить несколько раз с разными операндами, например, получить суммы нескольких пар чисел. SSE содержит команды для параллельного выполнения четырёх однотипных операций над вещественными числами одинарной точности. Появившийся чуть позже SSE2 содержит аналогичные команды для чисел двойной точности. Хотя SSE выполняет четыре операции одновременно, ничто не мешает использовать его там, где нужна только одна операция: результаты остальных трёх при этом просто игнорируются.
Первые версии .NET Framework были разработаны до широкого распространения SSE, так что эта платформа для операций с вещественными числами использует инструкции FPU, а появившийся позже .NET использует уже инструкции SSE. Этим и объясняется разница в результатах.
SSE выполняет вычисления «честно», не конвертируя свои аргументы в другой формат. Значение 0.1f загружается в соответствующий регистр как есть, а операция умножения на 10 даёт минимально возможную для данного типа погрешность результата, которая в данном случае удачно оказывается равной нулю.
FPU расширяет мантиссу числа одинарной точности до 63 знаков, заполняя их, как обычно, нулями. Дальнейшая работа идёт в пределах погрешности, которая достигается на 63 битах, т.е. существенно меньшей, чем погрешность одинарной точности. Как мы знаем, 0.1f это примерно 0.10000000149. Умножая на 10, получаем 1.0000000149, вычитая результат из 1, получаем 1.49e-8 – примерно такое число мы и видим на экране при использовании .NET Framework.
Итак, разница заключается в том, что .NET (точнее, SSE) работает с числами с той погрешностью, с которой они хранятся в памяти, а .NET Framework (точнее, FPU) использует промежуточное представление с существенно меньшей погрешностью, чем погрешность исходных данных. В случае .NET погрешность расчётов как бы компенсирует погрешность исходных данных, а в .NET погрешность расчёта намного ниже и не может компенсировать ту же самую погрешность исходных данных.
Но не надо думать, что для любого целого n в .NET будет выполняться равенство n*(1/n)=1. Это нам только с числом 0.1 так повезло, с другими числами полная компенсация погрешностей может и не произойти (хотя в большинстве случаев всё-таки произойдёт). Вот примеры этого.
Console.WriteLine("n*(1/n) != 1");
f = 1.0f / 47;
f = 1 - f * 47;
Console.WriteLine($"1 - 47 * (1 / 47) Float: {f}");
d = 1.0d / 47;
d = 1 - d * 47;
Console.WriteLine($"1 - 47 * (1 / 47) Double: {d}");
f = 1.0f / 49;
f = 1 - f * 49;
Console.WriteLine($"1 - 49 * (1 / 49) Float: {f}");
d = 1.0d / 49;
d = 1 - d * 49;
Console.WriteLine($"1 - 49 * (1 / 49) Double: {d}");
Результат в .NET:
n*(1/n) != 1
1 - 47 * (1 / 47) Float: 5.9604645E-08
1 - 47 * (1 / 47) Double: 0
1 - 49 * (1 / 49) Float: 0
1 - 49 * (1 / 49) Double: 1.1102230246251565E-16
Результат в .NET Framework:
n*(1/n) != 1
1 - 47 * (1 / 47) Float: 3.166497E-08
1 - 47 * (1 / 47) Double: 0
1 - 49 * (1 / 49) Float: 2.04891E-08
1 - 49 * (1 / 49) Double: 1.11022302462516E-16
Мы видим, что .NET Framework ожидаемо не справляется с числами одинарной точности, но и .NET справляется с ними тоже не всегда. А также существуют числа двойной точности, с которыми не справляется даже .NET.
Управляющее слово FPU
Если вы хорошо поняли, почему 10*0.1f в случае .NET Framework даёт не единицу, у вас должен появиться один вопрос: но ведь числа двойной точности тоже имеют погрешность, большую, чем погрешность внутреннего 10-байтного представления FPU, и там погрешность расчёта тоже не соответствует погрешности исходных данных. Почему же тогда мы получаем единицу? Это очень правильный вопрос, и чтобы ответить на него, нужно познакомиться ещё с некоторыми особенностями FPU.
У FPU есть специальный 16-битный регистр, называемый управляющим словом. Отдельные биты этого регистра позволяют задавать разные настройки, используемые при вычислениях. Среди прочих настроек там есть управление точностью вычислений. За это отвечают два бита – восьмой и девятый. 00 означает одинарную точность, 10 – двойную точность, 11 – расширенную (значение 01 зарезервировано).
По умолчанию FPU стартует в режиме расширенной точности, в котором используются все биты мантиссы. В режиме двойной точности используются только 53 бита мантиссы, что совпадает с размером мантиссы числа двойной точности IEEE 754 (строго говоря, у числа двойной точности 52 бита в мантиссе, но в формате FPU добавляется один бит для записи разряда перед точкой). В режиме одинарной точности используется, соответственно, 24 бита мантиссы. Я так и не нашёл внятного ответа на вопрос, зачем нужны такие ограничения точности, но, по некоторым сведениям, в старых версиях сопроцессоров это позволяло увеличить скорость вычислений.
Обратите внимание, что выбор соответствующего режима точности не означает, что FPU начинает работать с числами одинарной или двойной точности, так как, кроме мантиссы, есть ещё и экспонента, на размер которой данный режим не влияет.
.NET Framework переключает FPU в режим двойной точности. Я могу только догадываться, почему так сделано – возможно, для большей совместимости с аппаратными платформами, которые не поддерживают 10-байтный формат. Именно из-за этого результаты вычисления с числами двойной точности в .NET и .NET Framework обычно совпадают (дальше мы увидим, что могут и не совпадать).
.NET Framework, разумеется, не предоставляет возможности напрямую изменять управляющее слово FPU, но когда очень хочется, то можно. Сразу скажу, что способ, который я сейчас приведу, я крайне не рекомендую применять в реальных программах, но в исследовательских целях такое вполне допустимо. В состав Windows включены библиотеки времени выполнения для программ на C/C++, и мы можем использовать эти библиотеки. Нужная нам функция находится в библиотеке msvcrt.dll и называется _controlfp. Её можно импортировать следующим образом:
const int _MCW_PC = 0x00030000;
const int _PC_24 = 0x00020000;
const int _PC_53 = 0x00010000;
const int _PC_64 = 0x00000000;
[DllImport("msvcrt.dll", CallingConvention = CallingConvention.Cdecl)]
public static extern int _controlfp(int newControl, int mask);
Я не буду подробно описывать функцию _controlfp и её параметры, скажу только, что в нашем случае в качестве первого параметра надо передавать одну из констант _PC_xx (в зависимости от режима точности, который мы хотим выбрать), а в качестве второго параметра – константу _MCW_PC.
Обратите внимание, что код, приведённый в этом разделе, имеет смысл только для .NET Framework, так как в .NET изменение управляющего слова FPU не приведёт ни к какому результату. Соответственно, этот код есть только в примере для .NET Framework.
Попробуем повторить вычисления из предыдущего раздела, меняя режимы точности.
Console.WriteLine("0.1*10, 24 bits");
_controlfp(_PC_24, _MCW_PC);
f = 0.1f;
f = 1 - 10 * f;
Console.WriteLine($"1 - 0.1f*10: {f}");
d = 0.1;
d = 1 - 10 * d;
Console.WriteLine($"1 - 0.1d*10: {d}");
Console.WriteLine();
Console.WriteLine("0.1*10, 64 bits");
_controlfp(_PC_64, _MCW_PC);
f = 0.1f;
f = 1 - 10 * f;
Console.WriteLine($"1 - 0.1f*10: {f}");
d = 0.1;
d = 1 - 10 * d;
Console.WriteLine($"1 - 0.1d*10: {d}");
Результат в .NET Framework:
0.1*10, 24 bits
1 - 0.1f*10: 0
1 - 0.1d*10: 0
0.1*10, 64 bits
1 - 0.1f*10: -1.490116E-08
1 - 0.1d*10: -5.55111512312578E-17
При переводе FPU в режим 24-битной мантиссы при вычислениях с числами одинарной точности мы получаем 0, потому что теперь у нас нет разницы в точности исходных данных и расчёта, имеем тот же результат, что и в случае .NET. Что касается вычислений с двойной точностью, то тут при загрузке в регистр FPU число 0.1d округляется до 0.1f, с которым выполняются те же самые вычисления, что и в первом случае, поэтому получаем такой же результат. А вот в режиме 64-битной мантиссы мы даже для исходных данных двойной точности имеем погрешность вычислений, меньшую, чем погрешность исходных данных, что закономерно приводит к ненулевому результату.
Возможность изменения точности приводит к одному очень неприятному эффекту: если вы используете любой сторонний код, он может изменить точность вычислений по своему усмотрению, и ваши расчёты до и после вызова этого кода будут давать разные результаты. Особенно это касается внешних библиотек для сложных вычислений, разработчики которых могут иметь своё мнение насчёт того, какой режим им лучше подходит. Это маловероятно, если внешняя библиотека тоже написана для .NET Framework, но если это нативная библиотека, вероятность наткнуться на такую ситуацию возрастает. Если вы используете подобные библиотеки, рекомендую проверить, не меняют ли они управляющее слово FPU.
В дальнейшем в этой статье, если явно не оговорено иное, для .NET Framework приводятся результаты, соответствующие режиму 53-битной мантиссы, который устанавливается по умолчанию.
Вычитание в цикле
Выше мы умножали 0.1 на 10 и вычитали произведение из единицы. Проверим, изменится ли результат, если вместо умножения мы вычтем из 0.1 из единицы десять раз.
Console.WriteLine("Subtraction in cycle");
f = 1;
for (int i = 0; i < 10; i++)
f -= 0.1f;
Console.WriteLine($"1.0f - 0.1f (10 times): {f}");
d = 1;
for (int i = 0; i < 10; i++)
d -= 0.1;
Console.WriteLine($"1.0d - 0.1d (10 times): {d}");
Результат в .NET:
Subtraction in cycle
1.0f - 0.1f (10 times): -7.4505806E-08
1.0d - 0.1d (10 times): 1.3877787807814457E-16
Результат в .NET Framework:
Subtraction in cycle
1.0f - 0.1f (10 times): -7.450581E-08
1.0d - 0.1d (10 times): 1.38777878078145E-16
Хорошая новость: результаты .NET и .NET Framework совпали. Плохая новость: и там, и там они неправильные. Ожидаемого ноля мы не получили нигде.
Здесь мы видим результат накопления ошибки: при каждом действии имеем ошибку в одну сторону, в результате никакой компенсации не происходит.
А если вычитать не в цикле, а в одном выражении? Будет ли разница? Давайте проверим.
Console.WriteLine("Subtraction 10 times");
f2 = 0.1f;
f = 1 - f2 - f2 - f2 - f2 - f2 - f2 - f2 - f2 - f2 - f2;
Console.WriteLine($"1.0f - 0.1f (10 times): {f}");
d2 = 0.1;
d = 1 - d2 - d2 - d2 - d2 - d2 - d2 - d2 - d2 - d2 - d2;
Console.WriteLine($"1.0d - 0.1d (10 times): {d}");
Результат в .NET:
Subtraction 10 times
1.0f - 0.1f (10 times): -7.4505806E-08
1.0d - 0.1d (10 times): 1.3877787807814457E-16
Результат в .NET Framework:
Subtraction 10 times
1.0f - 0.1f (10 times): -1.490116E-08
1.0d - 0.1d (10 times): 1.38777878078145E-16
В .NET результаты не изменились, в .NET Framework изменился результат расчёта с числами одинарной точности. Чтобы понять причину этого, нужно выяснить, чем вычитание в одном выражении отличается от вычитания в цикле. Принципиальное отличие одно: при вычитании в цикле промежуточный результат каждого вычитания сохраняется в переменной, а при вычитании в одном выражении он сохраняется в регистрах процессора. В .NET используются регистры SSE, которые имеют ту же размерность, что и переменные, поэтому на результат это не влияет. А в .NET Framework используются 10-байтные регистры FPU, в которых промежуточный результат сохраняется с меньшей погрешностью, чем в переменной, и только окончательный результат усекается до размерности типа одинарной точности. Поэтому результат для вычислений с одинарной точностью получается другим.
Что касается вычислений с двойной точностью в .NET Framework, то напоминаю, что эта платформа устанавливает такой режим работы FPU, когда точность хранения значений в регистрах FPU совпадает с точностью типа double. Поэтому, как и в случае .NET, точность промежуточного результата не зависит от того, хранится он в переменной или в регистре. Но если вы переведёте FPU в режим 64-битной мантиссы, то результат вычислений с двойной точностью тоже изменится.
Но вернёмся к вычитанию в цикле. Результат, приведённый для .NET Framework, вы увидите только при компиляции проекта в конфигурации Debug. Если вы запустите проект в конфигурации Release, то для типа float при вычитании в цикле вы увидите то же значение -1.490116E-08, что и при вычитании без цикла. Это связано с тем, что в конфигурации Debug результирующий код делает в точности то, что написано в исходном тексте: если там стоит запись значения в переменную, то это значение реально будет занесено в соответствующие ячейки памяти, а при обращении к этой переменной – прочитаны из этих ячеек. А в конфигурации Release платформа старается, в том числе и при вещественных вычислениях, максимально уменьшить количество обращений к памяти и сохранять промежуточные результаты в регистрах процессора, так как это намного быстрее. В этой конфигурации на переменные, которые служат только для хранения промежуточных результатов, часто вообще не выделяется место в памяти, компилятор обходится регистрами. Соответственно, для вычитания в цикле промежуточный результат в конфигурации Release сохраняется не в переменной, а в регистре FPU. По сути, получается та же ситуация, что и при вычитании без цикла, и на выходе будет тот же результат.
В .NET конфигурация Release тоже старается обойтись регистрами вместо переменных везде, где это возможно. Но так как точность регистров SSE такая же, как у переменных, это не влияет на конечный результат.
Машинное эпсилон
В предыдущем разделе мы выяснили, что, вычтя 0.1 из единицы 10 раз, мы не получим 0. Начинающие программисты, которые этого не знают, иногда пытаются писать такой код:
float f = 1;
while (f != 0)
{
...
f -= 0.1f;
}
Каков будет результат выполнения такого цикла? Сразу хочется сказать, что вычитание будет продолжаться до тех пор, пока значение f не достигнет минус бесконечности, после чего перестанет меняться, так как при вычитании из минус бесконечности любого числа мы вновь будем получать минус бесконечность. Это неправильный ответ. На самом деле переменная f перестанет меняться, когда достигнет значения -2097152. Вычитая из этого числа 0.1, мы каждый раз будем снова получать -2097152. Это происходит потому, что одинарной точности не хватает, чтобы записать число -2097152.1: следующее возможное число – это -2097152.2. Результат округляется вверх, и мы получаем -2097152, с которого начали. Другими словами, добавляемое число меньше погрешности, поэтому результат равен первому слагаемому.
Для характеристики этой особенности вычислений с погрешностью вводится понятие машинного эпсилон: это минимальное положительное число, при добавлении которого к единице получается значение, отличное от единицы. При добавлении к единице любого числа, которое по модулю меньше машинного эпсилон, мы получим единицу.
Давайте вычислим, чему равно машинное эпсилон. Единица это 1*2^0, т.е. для его представления мантисса должна содержать единицу, а экспонента – значение, соответствующее нулевой степени. Так как целая единица в нормализованной записи только подразумевается, это значит, что все биты мантиссы будут равны нулю. Соответственно, минимальное отличное от единицы число мы получим, если при той же экспоненте сделаем единицей самый младший бит мантиссы. Машинное эпсилон будет разностью этого числа и единицы. Таким образом, для чисел одинарной точности машинное эпсилон будет равно 2^-23 = 1.19e-7, для чисел двойной точности – 2^-52=2.22e-16, а для 10-байтных регистров FPU – 2^-63=1.08e-19. Проверим наши расчёты практикой.
Console.WriteLine("Simple epsilon search");
f = 1;
while (1 + f / 2 > 1)
f /= 2;
Console.WriteLine($"Float epsilon: {f}");
d = 1;
while (1 + d / 2 > 1)
d /= 2;
Console.WriteLine($"Double epsilon: {d}");
Результат в .NET:
Simple epsilon search
Float epsilon: 1.1920929E-07
Double epsilon: 2.220446049250313E-16
Результат в .NET Framework:
Simple epsilon search
Float epsilon: 2.220446E-16
Double epsilon: 2.22044604925031E-16
Результаты в .NET соответствуют нашим ожиданиям, а в .NET Framework значения для типов float и double оказались равны (они действительно в точности равны друг другу, более точное значение для double – это особенности работы метода ToString, который выводит разные количества знаков для float и double).
Вам уже вполне хватает знаний, чтобы самостоятельно догадаться, почему так получилось. Мы сталкивались с аналогичной ситуацией в предыдущем разделе. Промежуточные результаты вычисления выражения 1+f/2>1 хранятся в регистрах процессора, и в случае .NET Framework это регистры FPU, имеющие в режиме по умолчанию точность типа double. Поэтому вычисления производятся с точностью этого типа и дают в итоге тот же результат. Избежать этого можно, если сохранять промежуточные результаты в переменной требуемого типа.
Console.WriteLine("Epsilon search by steps");
f = 1;
while (true)
{
f2 = f / 2;
f2 = 1 + f2;
if (f2 == 1)
break;
f /= 2;
}
Console.WriteLine($"Float epsilon: {f}");
d = 1;
while (true)
{
d2 = d / 2;
d2 = 1 + d2;
if (d2 == 1)
break;
d /= 2;
}
Console.WriteLine($"Double epsilon: {d}");
Результат в .NET:
Epsilon search by steps
Float epsilon: 1.1920929E-07
Double epsilon: 2.220446049250313E-16
Результат в .NET Framework:
Epsilon search by steps
Float epsilon: 1.192093E-07
Double epsilon: 2.22044604925031E-16
Теперь в обеих платформах получаем правильный результат.
Замечу, что в .NET Framework в конфигурации Release компилятор в этом коде опять обойдётся регистрами FPU, поэтому даже код с промежуточной переменной выдаст для float то же значение, что и для double. Правильное значение для float удаётся получить только в конфигурации Debug.
Если при простом вычислении машинного эпсилон в .NET Framework результат зависит не от точности исходных данных, а от точности регистров FPU, изменение этой точности должно влиять на результат. Давайте это проверим. Я не буду приводить здесь код для вычисления, потому что он полностью аналогичен тому, что уже есть выше, только скомпонован немного по-другому, и добавлено переключение размерности мантиссы. Код есть в примере FPNumbersFW.
Результат в .NET Framework:
Epsilon, 24 bits
Float epsilon, simple: 1.192093E-07
Float epsilon, by steps: 1.192093E-07
Double epsilon, simple: 1.19209289550781E-07
Double epsilon, by steps: 1.19209289550781E-07
Epsilon, 64 bits
Float epsilon, simple: 1.084202E-19
Float epsilon, by steps: 1.192093E-07
Double epsilon, simple: 1.0842021724855E-19
Double epsilon, by steps: 2.22044604925031E-16
В режиме 24-битной мантиссы, как ни вычисляй эпсилон и какие типы при этом ни используй, ограничивающим фактором является точность регистров, поэтому мы везде получаем значение эпсилон, соответствующее одинарной точности. В режиме 64-битной мантиссы при простом вычислении точность ограничена точностью регистров, и мы получаем значение, соответствующее максимальной точности регистров. При вычислении по шагам ограничивающим фактором становится точность переменной, в которую записывается промежуточное значение, поэтому мы получаем значения, соответствующие точности выбранных переменных.
Перестановка слагаемых
При перестановке слагаемых сумма не изменяется – это правило мы помним с первого класса. Но для вычислений с погрешностью оно не всегда работает. Когда речь идёт о двух слагаемых, то их можно складывать в любом порядке, но если слагаемых больше, то тут уже возможны варианты, и чем больше слагаемых, тем больше вероятность попасть в такую ситуацию. Рассмотрим следующий код.
Console.WriteLine("Terms rearrangement");
double[] da = new double[1000000];
for (int i = 1; i < da.Length; i++)
da[i] = 100;
da[0] = 1e19;
d = 0;
for (int i = 0; i < da.Length; i++)
d += da[i];
Console.WriteLine($"Direct sum: {d}");
d = 0;
for (int i = da.Length - 1; i >= 0; i--)
d += da[i];
Console.WriteLine($"Reverse sum: {d}");
Результат в .NET:
Terms rearrangement
Direct sum: 1E+19
Reverse sum: 1.00000000001E+19
Результат в .NET Framework:
Terms rearrangement
Direct sum: 1E+19
Reverse sum: 1.00000000001E+19
Возможно, вы уже догадались, откуда берётся разница. Число 100 меньше, чем погрешность записи числа 1e19, поэтому, когда мы прибавляем 100 к 1e19, то получаем снова 1e19, и от того, что мы повторили эту операцию почти миллион раз, результат не изменится. А вот когда мы сначала складываем почти миллион раз число 100 с самим собой, получившийся результат уже больше погрешности 1e19, поэтому при сложении мы получаем уже не 1e19.
Ещё более интересный эффект можно получить, если заменить в этом коде 1e19 на 1e18.
Результат в .NET:
Terms rearrangement (with 1e18)
Direct sum: 1.0000000001279999E+18
Reverse sum: 1.0000000000999999E+18
Результат в .NET Framework:
Terms rearrangement (with 1e18)
Direct sum: 1.000000000128E+18
Reverse sum: 1.0000000001E+18
Опять замечу, что результаты на двух платформах одинаковые, просто в .NET Framework метод ToString округляет сильнее.
Число 100 превышает погрешность числа 1e18, поэтому, когда мы их складываем, получается уже не 1e18, а другое число. Но число 100 в двоичном представлении не круглое, это 64+32+4, т.е. 2^6+2^5+2^2. При записи этого числа с плавающей двоичной точкой в экспоненте будет значение, соответствующее шести, а мантисса будет выглядеть как [1].100100…, где подразумеваемый целый разряд даёт (с учётом значения экспоненты) 64, первый дробный разряд – 32 и т.д.
Теперь представьте, как 100 прибавляется к 1e18. Тут так получилось, что самый старший, подразумеваемый двоичный разряд числа 100 по порядку соответствует самому младшему разряду мантиссы числа 1e18. Всё, что идёт в числе 100 после двоичной точки, меньше погрешности представления 1e18, так что прибавить эти разряды к 1e18 невозможно. Но они не просто отбрасываются, а округляются вверх, так как 0.1001 – это больше, чем половина. Соответственно, получается, что на каждом шаге к 1e18 прибавляется не 100, а 128, что мы и видим в значении результата.
Проверить это утверждение можно очень простым кодом:
Console.WriteLine("1e18 and 100");
d = 1e18 + 100 - 1e18;
Console.WriteLine($"1e18 + 100 - 1e18: {d}");
Результат в .NET:
1e18 and 100
1e18 + 100 - 1e18: 128
Результат в .NET Framework:
1e18 and 100
1e18 + 100 - 1e18: 128
Тут есть один интересный момент: на самом деле компилятор не создаёт код, который складывает 1e18 и 100, а потом вычитает 1e18. Это константное выражение, которое вычисляется на этапе компиляции, а в коде в переменную d сразу записывается результат вычисления. Но так как для вычисления выражения компилятор пользуется той же вещественной арифметикой, то и результат у него получается такой же, как если бы выражение вычислялось в коде программы.
Перестановка множителей
Результат умножения тоже может зависеть от порядка множителей. Вот пример:
Console.WriteLine("Factors rearrangement");
d = 1e200;
d = d * d * 1e-300;
Console.WriteLine($"1e200*1e200*1e-300: {d}");
d = 1e200;
d = d * 1e-300 * d;
Console.WriteLine($"1e200*1e-300*1e200: {d}");
d = 1e-200;
d = d * d * 1e300;
Console.WriteLine($"1e-200*1e-200*1e300: {d}");
d = 1e-200;
d = d * 1e300 * d;
Console.WriteLine($"1e-200*1e300*1e-200: {d}");
Результат в .NET:
Factors rearrangement
1e200*1e200*1e-300: ∞
1e200*1e-300*1e200: 1E+100
1e-200*1e-200*1e300: 0
1e-200*1e300*1e-200: 1E-100
Результат в .NET Framework:
Factors rearrangement
1e200*1e200*1e-300: 1E+100
1e200*1e-300*1e200: 1E+100
1e-200*1e-200*1e300: 1E-100
1e-200*1e300*1e-200: 1E-100
Мы видим, что в .NET результат зависит от порядка выполнения операций, а в .NET Framework не зависит. На самом деле, и в .NET Framework результат в некоторых случаях может зависеть от порядка множителей, но я специально подобрал такие значения, чтобы продемонстрировать ещё один интересный эффект, из-за которого результаты вычислений в .NET и .NET Framework могут различаться.
Сначала разберёмся, что происходит в .NET. В результате умножения 1e200 на 1e200 должно получиться 1e400, но это значение выходит за пределы двойной точности (там максимум, напомню, 1.8e308), поэтому оно заменяется на бесконечность. Когда мы умножаем бесконечность на конечное число 1e-300, снова получаем бесконечность. А вот если переставить множители, то промежуточный вариант не выходит за пределы допустимого диапазона, и мы получаем правильное значение. Аналогично при умножении 1e-200 на 1e-200 мы получаем число 1e-400, которое слишком мало для чисел двойной точности и потому заменяется на ноль. При перестановке множителей выхода за пределы двойной точности не происходит.
В .NET Framework всё дело, как обычно, в 10-байтных регистрах, в которых сохраняются промежуточные результаты. Диапазон этих регистров – от 3.6e-4951 до 1.18e4932, поэтому и 1e400, и 1e-400 помещаются туда без проблем. Соответственно, эффекта от перестановки множителей мы не видим. При использовании FPU, конечно, тоже можно добиться выхода за пределы диапазона, но сделать это в .NET Framework не очень просто, так как исходные данные могут иметь максимум двойную точность. Даже максимальное значение двойной точности, возведённое в квадрат, даёт всего лишь 3.2e616, что с большим запасом помещается в диапазон регистров FPU. Требуется многократное перемножение чисел двойной точности, чтобы выйти за пределы диапазона 10-байтных регистров.
Ещё обратите внимание, что результат работы примера в .NET Framework не зависит от того, какой режим мантиссы выбран. Этот режим накладывает ограничения только на размер мантиссы, не затрагивая экспоненту, а за диапазон отвечает именно экспонента.
Округление
За округление к ближайшему значению в .NET (Framework) отвечает метод Math.Round. Однако, если вы никогда не писали приложений для финансовой сферы и не знакомы с принятыми в этой области правилами округления, ваши ожидания по поводу результатов работы этой функции, скорее всего, неверны. Вот простой код, который демонстрирует эти особенности.
Console.WriteLine("Rounding");
Console.WriteLine($"[2.5] - [1.5]: {Math.Round(2.5) - Math.Round(1.5)}");
Console.WriteLine($"[3.5] - [2.5]: {Math.Round(3.5) - Math.Round(2.5)}");
Результат в .NET:
Rounding
[2.5] - [1.5]: 0
[3.5] - [2.5]: 2
Результат в .NET Framework:
Rounding
[2.5] - [1.5]: 0
[3.5] - [2.5]: 2
Если бы метод Math.Round выполнял обычное арифметическое округление, то в обоих случаях мы бы увидели единицу. Но данный метод использует другой вариант округления, известный как банковское, или бухгалтерское, округление.
Вообще, стандарт IEEE 754 описывает пять видов округления, и это даже не полный список того, что можно придумать. Округление к нулю – из двух чисел, между которыми находится округляемое число, выбирается то, которое меньше по модулю, т.е. ближе к нулю (данный вид округления в .NET (Framework) реализуется методом Math.Truncate). Округление к плюс бесконечности – выбирается то число, которое больше (метод Math.Ceiling). Округление к минус бесконечности – выбирается то число, которое меньше (метод Math.Floor). Округление к ближайшему с привязкой к бесконечности – из двух чисел выбирается то, которое ближе к округляемому, а если оно находится точно между ними, то выбирается то, которое больше по модулю (это и есть обычное арифметическое округление). Округление к ближайшему с привязкой к чётному – если округляемое значение лежит точно между двумя значениями, из них выбирается чётное значение. Последний вариант и называется банковским округлением, и именно его реализует метод Math.Round с одним параметром. Поэтому и 1.5, и 2.5 округляются этим методом до 2, а 3.5 – до 4, что приводит к тем результатам, которые мы видим.
Пара слов о том, зачем это понадобилось банкам и бухгалтерам. Им часто приходится решать задачи вида «три человека выполнили работу общей стоимостью 100 рублей, внеся равные вклады, как разделить заработанные деньги между ними поровну?» Конечно, математику не обманешь и 100 рублей на троих поровну не разделишь, но у бухгалтеров есть свои методы решения подобных проблем (не спрашивайте меня, какие именно, в этом я плохо разбираюсь). И эти методы плохо стыкуются с арифметическим округлением, потому что для него математическое ожидание суммы округлённых величин не равно сумме величин до округления. Грубо говоря, ошибка, вносимая округлением значения 0.1, компенсируется такой же по модулю, но противоположной по знаку ошибкой, возникающей при округлении 0.9, 0.2 компенсируется 0.8 и т.д., а вот 0.5 ничем не компенсируется, всегда округляется вверх по модулю. Банковское округление лишено этого недостатка – там чётное 0.5 компенсируется нечётным 0.5, и бухгалтерам такое округление нравится намного больше.
Существует перегрузка метода Math.Round, которая позволяет явно выбрать любой из стандартных способов округления, в т.ч. арифметическое округление:
Console.WriteLine("Arithmetic rounding");
Console.WriteLine($"[2.5] - [1.5]: {Math.Round(2.5, MidpointRounding.AwayFromZero) - Math.Round(1.5, MidpointRounding.AwayFromZero)}");
Console.WriteLine($"[3.5] - [2.5]: {Math.Round(3.5, MidpointRounding.AwayFromZero) - Math.Round(2.5, MidpointRounding.AwayFromZero)}");
Результат в .NET:
Arithmetic rounding
[2.5] - [1.5]: 1
[3.5] - [2.5]: 1
Результат в .NET Framework:
Arithmetic rounding
[2.5] - [1.5]: 1
[3.5] - [2.5]: 1
Ещё один момент, о котором часто забывают, – это то, что методы Math.Ceiling и Math.Floor несимметричны относительно нуля, т.е. для любого нецелого x верно неравенство Math.Ceiling(x) != -Math.Ceiling(-x), а также аналогичное неравенство для Math.Floor. Следующий пример это показывает.
Console.WriteLine("Floor and ceiling");
Console.WriteLine($"Floor [2.5]: {Math.Floor(2.5)}");
Console.WriteLine($"Floor [-2.5]: {Math.Floor(-2.5)}");
Console.WriteLine($"Ceiling [2.5]: {Math.Ceiling(2.5)}");
Console.WriteLine($"Ceiling [-2.5]: {Math.Ceiling(-2.5)}");
Результат в .NET:
Floor and ceiling
Floor [2.5]: 2
Floor [-2.5]: -3
Ceiling [2.5]: 3
Ceiling [-2.5]: -2
Результат в .NET Framework:
Floor and ceiling
Floor [2.5]: 2
Floor [-2.5]: -3
Ceiling [2.5]: 3
Ceiling [-2.5]: -2
У метода Math.Round существуют также перегрузки, которые позволяют округлять не до целого, а до какого-то десятичного знака. Но я бы рекомендовал относиться к этой возможности скептически. Во-первых, как мы уже знаем, большинство десятичных дробей представляется числами с плавающей двоичной точкой только приближённо, поэтому округление до этих чисел – это само по себе уже некоторая условность. Во-вторых, само округляемое число тоже представлено с некоторой ошибкой, которая может быть как положительной, так и отрицательной. Вполне возможна ситуация, что число, которое оканчивается на 5, в машинном представлении чуть меньше, чем в десятичной записи. Тогда арифметическое округление округлит его вниз, а не вверх. Это иллюстрируется следующим примером.
Console.WriteLine("Rounding to hundredths");
for (int i = 0; i <= 9; i++)
{
d = double.Parse($"2.1{i}5");
Console.WriteLine($"{d}: {Math.Round(d, 2, MidpointRounding.AwayFromZero)}");
}
Результат в .NET:
Rounding to hundredths
2.105: 2.11
2.115: 2.12
2.125: 2.13
2.135: 2.13
2.145: 2.15
2.155: 2.15
2.165: 2.17
2.175: 2.17
2.185: 2.19
2.195: 2.19
Результат в .NET Framework:
Rounding to hundredths
2.105: 2.11
2.115: 2.12
2.125: 2.13
2.135: 2.13
2.145: 2.15
2.155: 2.15
2.165: 2.17
2.175: 2.17
2.185: 2.19
2.195: 2.19
Во взятом наугад диапазоне из десяти чисел четыре округлены неправильно.
Ещё более интересные эффекты можно получить при округлении к плюс бесконечности или к минус бесконечности. Число, которое уже как бы округлено, может превратиться не само в себя, а немного измениться. Такие эффекты достаточно редки, потому что операция округления, в целом, учитывает эти тонкости, но всё-таки они встречаются.
Console.WriteLine("Ceiling rounding to hundredths");
for (int i = 0; i <= 9; i++)
{
d = double.Parse($"2.1{i}");
Console.WriteLine($"{d}: {Math.Round(d, 2, MidpointRounding.ToPositiveInfinity)}");
}
Результат в .NET:
Ceiling rounding to hundredths
2.1: 2.1
2.11: 2.11
2.12: 2.12
2.13: 2.13
2.14: 2.14
2.15: 2.15
2.16: 2.16
2.17: 2.17
2.18: 2.19
2.19: 2.19
Результат приводится только для .NET, потому что в .NET Framework этот код невоспроизводим: метод Math.Round там реализует только арифметическое и банковское округление, а метод Math.Ceiling умеет округлять только до целых.
Здесь мы неверный результат для числа 2.18, которое превратилось в 2.19, остальные числа остались без изменений.
Тип decimal
В .NET (Framework) есть ещё один тип для хранения вещественных чисел – System.Decimal (decimal в C#). Я не буду здесь подробно останавливаться на этом типе, так как он несколько выходит за рамки этой статьи, но некоторые общие сведения приведу, чтобы дать хотя бы общее представление о том, чем он принципиально отличается от float и double.
Значение типа decimal занимает в памяти 8 байт – 128 бит. Из них 1 бит отводится на знак, 96 – на целочисленное базовое значение B, 8 бит отводится на степень делителя F, 23 бита не используются. При этом степень делителя должна находится в диапазоне от 0 до 28, хотя 8 бит могут хранить значение до 255. Модуль значения вычисляется по формуле B/10^F. Т.е., например, для записи значения 12345 B=12345, F=0. Значение 123.45 – это B=12345 и F=2. А 1.2345 – это B=12345 и F=4.
Описанный выше способ кодирования числа избыточен, т.е. одно число может быть записано несколькими способами. Например, 123.45 это ещё и B=123450, F=3, а также B=1234500, F=4 и т.д. Все эти способы записи одного числа имеют право на существование. Это затрудняет реализацию операции сравнения, потому что она должна учитывать равенство всех подобных вариантов. Впрочем, эти сложности скрыты внутри платформы, для программиста, использующего тип decimal, это абсолютно прозрачно.
Так как decimal работает с десятичными, а не двоичными степенями, у него нет проблем с точным представлением десятичных дробей до 28 знаков. Скажем, с задачей на вычитание 0.1 из 1 десять раз, на которой спотыкаются числа с плавающей двоичной точкой, тип decimal справляется без проблем.
Console.WriteLine("Decimal subtraction");
decimal m = 1;
decimal m2 = 0.1m;
for (int i = 0; i < 10; i++)
m -= m2;
Console.WriteLine($"1 - 0.1m (10 times): {m}");
Результат в .NET:
Decimal subtraction
1 - 0.1m (10 times): 0.0
Результат в .NET Framework:
Decimal subtraction
1 - 0.1m (10 times): 0.0
Платой за точность является более медленная скорость работы. Тип decimal не имеет аппаратной реализации, операции с ним выполняются подпрограммами платформы, поэтому работа с ним существенно медленнее. Чтобы понять, насколько, сделаем замеры (в многозадачной системе эти замеры, понятное дело, весьма неточны, но для грубой оценки этого хватит).
Console.WriteLine("Performance measurement");
Stopwatch sw = Stopwatch.StartNew();
for (int i = 0; i < 10000000; i++)
{
d = 0.5;
d *= 10;
}
sw.Stop();
Console.WriteLine($"Time for double: {sw.ElapsedMilliseconds}");
sw = Stopwatch.StartNew();
for (int i = 0; i < 10000000; i++)
{
m = 0.5m;
m *= 10;
}
sw.Stop();
Console.WriteLine($"Time for decimal: {sw.ElapsedMilliseconds}");
Результат в .NET:
Performance measurement
Time for double: 10
Time for decimal: 225
Результат в .NET Framework:
Performance measurement
Time for double: 10
Time for decimal: 274
Результаты немного гуляют от запуска к запуску, но общее отношение примерно сохраняется (при запуске на другом компьютере, числа могут получиться совсем другими). Операция с типом decimal в двадцать с лишним раз медленнее, чем с типом double. При общей высокой производительности современных компьютеров для не слишком сложных вычислительных задач это не так уж и страшно, сотня-другая операций всё равно выполнятся за такое время, что пользователь не увидит разницы. Но для каких-то сложных расчётов типа математического моделирования разница может оказаться существенной.
На моём компьютере разница в скорости между .NET и .NET Framework в случае работы с double на уровне погрешности, а вот в случае decimal .NET стабильно показывает несколько лучший результат. Видимо, при разработке .NET подпрограммы для типа decimal были оптимизированы.
Если вы посмотрите стандарт IEEE 754, то увидите там типы с названием decimal, в частности, decimal128. Не надо путать его с типом System.Decimal. Кроме названия и размера, у этих типов нет ничего общего. Типы decimal в IEEE 754 тоже придуманы для хранения десятичных дробей без погрешностей, но решают эту задачу они совсем по-другому. Это тоже формат с плавающей точкой, но десятичной: мантисса там записывается с помощью двоично-десятичного кодирования, а экспонента хранит степени десяти, а не двух. Я не буду вдаваться в тонкости различных вариантов двоично-десятичных кодов (при том что IEEE 754 допускает использование двух вариантов этого кодирования), расскажу только самую простую идею. Для хранения числа от 0 до 9 требуется не менее четырёх бит, при этом из 16 возможных комбинаций этих бит 6 остаются невостребованными. Чтобы записать двузначное десятичное число, требуется 8 бит – по 4 бита на каждый разряд, трёхзначное – 12 бит и т.п. Это я описал самый простой способ двоично-десятичного кодирования, известный как BCD. В IEEE 754 применяются более продвинутые способы, которые позволяют записывать десятичные числа компактнее, но смысл остаётся тем же. Так как в мантиссе при этом, фактически, хранятся десятичные разряды, проблема с десятичным представлением дробей решена.
Полноценной поддержки десятичных форматов IEEE 754 в процессорах x86 и x64 нет. FPU формально их поддерживает, но только в смысле загрузке таких чисел из памяти в регистры и выгрузке их из регистров в память. При этом происходит преобразование десятичных форматов в уже знакомый нам 10-байтный формат с плавающей двоичной точкой, и точность представления десятичных дробей при этом теряется, т.е. точного расчёта с десятичными дробями FPU сделать не может. Эта возможность FPU чаще всего используется для преобразования вещественного числа в строку и наоборот, потому что десятично-двоичное число гораздо проще преобразовать в строку в виде десятичной дроби.
.NET (Framework) не имеет типов, соответствующих десятичным типам IEEE 754.
Заключение
Мы познакомились с основными ситуациями, когда вычисления с вещественными числами дают неожиданный результат. Конечно, большинство этих ситуаций связаны с тем, что вещественные числа используются на пределе возможностей, что редко бывает в реальной жизни. Тем не менее, шансы столкнуться с чем-то подобным всё-таки есть, и лучше заранее знать, что может случиться и почему.
Следует заметить, что в .NET вычисления ведут себя более понятно и предсказуемо, чем в .NET Framework, потому что в .NET Framework к проблемам, обусловленным особенностями стандарта IEEE 754, добавляются проблемы, причиной которых является внутренний формат представления чисел FPU и преобразование этого формата в форматы IEEE 754 и обратно. Особенно неприятной проблемой является то, что разные реализации одного и того же алгоритма могут давать разные результаты в зависимости от того, используются ли переменные для хранения промежуточного результата или нет. А с учётом склонности компилятора в конфигурации Release заменять переменные регистрами мы приходим к потенциально опасной ситуации, когда программа в режимах Debug и Release даёт разные результаты.
Я думаю, что отказ от FPU в пользу SSE разработчики .NET сделали именно потому, что хотели избежать проблем, связанных исключительно с FPU, и сделать код переносимым на другие платформы. Действительно, те результаты, которые обусловлены особенностями FPU, не могут быть получены на других аппаратных платформах, в то время как SSE полностью следует стандарту, и все результаты расчётов должны воспроизводиться на любой платформе, которая так же строго следует стандарту.
Перевод FPU в режим 53-разрядной мантиссы в .NET Framework, как я понимаю, тоже вызван желанием свести к минимуму те результаты, которые не могут быть получены в рамках чистого IEEE 754. Это имеет некоторый смысл, мы видели, что в этом режиме результаты в .NET Framework больше всего похожи на результаты в .NET. Но, во-первых, всё это касается только типа double, расчётам с использованием float это не помогает, а во-вторых, даже с double удаётся убрать далеко не все расхождения.
Избавляться от использования FPU разработчики платформы начали ещё до появления .NET. Если пример для .NET Framework откомпилировать в конфигурации x64, то он даст почти такие же результаты, как и пример для .NET. Отличия будут только в том, что 1^NaN даст NaN, а не 1, метод ToString для вещественных чисел будет выдавать меньше разрядов, чем в .NET, и сохранится более медленная работа с типом decimal. В 64-разрядном варианте .NET Framework тоже используется SSE вместо FPU. С одной стороны, это плохая новость: ваше приложение может начать работать иначе при переходе от 32-разрядной версии к 64-разрядной. С другой стороны, хорошая: если у вас 64-разрядное приложение, то его перенос в .NET, скорее всего, не доставит вам проблем в части вещественных вычислений.
Хуже обстоит дело с конфигурацией Any CPU. Согласно документации, приложения, откомпилированные в этой конфигурации, должны запускаться как 32-разрядные в 32-разрядных системах и как 64-разрядные в 64-разрядных системах. Однако в моей вполне 64-разрядной системе пример FPNumbersFW, будучи откомпилированным для Any CPU, запускается как 32-разрядное приложение и использует FPU. Предполагаю, что это может быть связано с маленьким размером приложения (.NET, кстати, ведёт себя иначе: FPNumbers в конфигурации Any CPU запускается как 64-разрядное приложение). Если моё предположение верно, то тут скрывается ещё один подводный камень: вы начинаете разрабатывать приложение, первое время оно невелико и запускается как 32-разрядное, но начиная с какого-то момента дорастает до 64-разрядного. И в этот момент у вас может перестать работать код, который вы считали уже отлаженным, а причина этого будет вам не очевидна.
Update: в .NET Framework для конфигурации Any CPU в настройках проекта есть опция "Предпочтительна 32-разрядная версия", которая по умолчанию включена. Видимо, из-за неё проект, откомпилированный для Any CPU, запускается в 32-разрядном режиме.
.NET и в этом отношении ведёт себя более предсказуемо: мне неизвестно ни одно случая, когда изменение разрядности или переключение между Debug и Release могли бы повлиять на результат вещественных операций. Более того, из всего, что я об этом знаю, следует, что таких ситуаций и не должно быть.
Но платформа .NET имеет реализации не только для x86 (64) и Windows, и то, как ведут себя эти реализации в других системах, требует отдельного исследования. Я таких исследований не проводил, так как работаю только с Windows, единственное, на что меня хватило, это запустить пример FPNumbers под WSL с установленной там Ubuntu. Я заметил только одно несущественное отличие: метод ToString место символа «∞» пишет слово «Infinity». Но вообще, вопрос воспроизводимости результатов при переносе на другие системы и аппаратные платформы очень интересный, если вы можете протестировать FPNumbers где-то ещё, напишите, пожалуйста, о результатах тестирования в комментариях.
Комментарии (21)
Krey
14.09.2023 11:23+1Монументально, большое спасибо.
Могу добавить что есть ещё библиотечный тип half
adeshere
14.09.2023 11:23Спасибо за фундаментальный обзор! Но у меня есть глупый вопрос:
> .NET Framework, разумеется, не предоставляет возможности напрямую изменять управляющее слово FPU
А как это слово инициализируется при старте программы? Разве ключи компиляции не позволяют это настроить? И если вдруг да, то как это на практике может работать? Например, что произойдет, если вдруг сторонняя функция изменит предустановленные настройки уже в рантайме?
Оффтопик с соседней улицы
У меня в устаревшем компиляторе Интел (фортран, 2013г) есть вот такой набор настроек FP-вычислений.
Сказать по правде, я никогда особо не задумывался над тем, что они означают (жил по принципу
работает- не трогай
Впрочем, один раз потрогать все же пришлось, когда оптимизатор неожиданно накосячил и
"забыл" вынуть из FPU результат вычисления
Я в одном месте вызывал функцию, которая возвращает результат FP-вычисления, но потом этот результат не использовал (от вызова функции мне нужен был побочный эффект). Пока я все собирал в режиме debug (без оптимизаций), все работало. А при включенной оптимизации начинало глючить - но очень редко и странно. Например генератор случайных чисел примерно один раз в 10 млн. вызовов вместо числа возвращал Nan.
Как оказалось, оптимизатор соображал, что результат вычисления не нужен, и доставать его незачем, но стек FPU при этом не чистил, запутавшись в море заинлайненых им же функций и их фрагментов. Отлаживать очень весело было, так как баг возникал вовсе не в момент вызова этой функции, а
совсем в другом месте при следующем обращении к FPU
Даже не знаю, сколько я слез пролил с этим багом, пока один мудрый человек не посоветовал включить контроль FPU-стека. После чего мое "UB" локализовалось в нужной точке, и дальше все стало просто.
Но после чтения Вашей статьи возникло подозрение, что на некоторые из описанных в ней эффектов (о которых я в общем-то знал, но избегал неприятностей другими способами) можно повлиять с помощью этих настроек?
Или нет?
Vglk Автор
14.09.2023 11:23+1С фортраном никогда не имел дела, но по моему небольшому опыту работы с интеловскими компиляторами для C++ я знаю, что они нацелены на то, чтобы выжать максимум именно из интеловских процессоров. В то время как компиляторы Microsoft стараются не использовать функции, доступные только на этих процессорах, чтобы код был более переносимый. Как пример, в C++ есть тип long double, который, согласно спецификации языка, совпадает с double или обеспечивает более высокую точность. В компиляторах Microsoft long double - это double, в компиляторах Intel это 80-битный тип, соответствующий по формату регистрам FPU.
> А как это слово инициализируется при старте программы?
У .NET-библиотек, которые отображаются в адресное пространство компилятора, есть секции инициализации, которые выполняются в момент этого отображения. В одной из них прописан перевод FPU в нужный режим.
> Разве ключи компиляции не позволяют это настроить?
У C# традиционно очень мало настроек компилятора. Ничего, связанного с вещественными числами, я в этих настройках не обнаружил.
> Например, что произойдет, если вдруг сторонняя функция изменит предустановленные настройки уже в рантайме?По сути, функция _controlfp, использующаяся в примерах, это и есть тот случай, когда сторонняя функция изменяет настройки. Возможно, они вернутся обратно после вызова какого-то метода .NET. По опыту работы с нативными приложениями знаю, что некоторые функции WindowsAPI (например, CreateWindow) переводят FPU в режим 53-битной мантиссы, поэтому результаты расчётов до вызова этой функции и после могут различаться.
adeshere
14.09.2023 11:23+1Как пример, в C++ есть тип long double, который, согласно спецификации языка, совпадает с double или обеспечивает более высокую точность. В компиляторах Microsoft long double - это double, в компиляторах Intel это 80-битный тип, соответствующий по формату регистрам FPU.
Боюсь, я совсем от жизни отстал. Я по наивности думал, что если я у себя в фортране напишу что-то вроде:
real (kind = 16) :: my_internal_real_var
real (kind = c_long_double) :: real_var_to_passто это как правило будет одно и то же, то есть 128-битный real. Даже прикидывал, как я этим смогу воспользоваться когда-нибудь в будущем...Получается, что все ровно наоборот, и как правило это будет совсем не одно и то же?
Расставаться с иллюзиями (что "мир устроен правильно и мудро" (с)) всегда неприятно, но все равно спасибо, что вернули мечтателя к суровой реальности ;-)
UPD:
Ну да, все печально - страница из справки фортрана
т.е. наличие прямых соответствий не только от компилятора зависит, но еще и от ОС...Просто от нас, программистов на фортране,
все это скрывают(пока не понадобится связка с другим языком)...Vglk Автор
14.09.2023 11:23Вообще, простые типы делятся на фундаментальные (base types) и обобщённые (generic types; не путать с generics - обобщениями, т.е. параметризованными типами). Фундаментальные типы - это типы, которые в любой реализации остаются одинаковыми, а обобщённые - те, которые могут меняться от реализации к реализации, чтобы наилучшим образом соответствовать конкретной среде. Так, в Паскале типы SmallInt и LongInt фундаментальные, это всегда и везде будет 16-битное и 32-битное целые соответственно. А тип Integer - обобщённый, в 16-разрядных версиях он был 16-битным, в 32-разрядных стал 32-битным. Не знаю, как с этим обстоят дела конкретно в Фортране, но общий подход сейчас такой.
adeshere
14.09.2023 11:23простые типы делятся на фундаментальные и обобщённые (...). Не знаю, как с этим обстоят дела конкретно в Фортране, но общий подход сейчас такой.
В старом фортране часто писали просто "integer", и это было похоже на паскалевский Integer (если не указать при компиляции разрядность такого "обобщенного" типа явно, то умолчание могло зависеть от реализации). Сейчас такой стиль тоже поддерживается для совместимости, но чаще все-таки пишут количество байт прямым текстом, например "INTEGER(KIND=8)". В теории, это должно обеспечить лучшую переносимость программ.
Кстати, 128-битных целых (KIND=16) в интеловском фортране 2013 года нет, только real.
a-tk
14.09.2023 11:23У компилятора C# нету опций относительно срезания углов с вещественными числами, потому что он следует букве стандарта IEEE754. А компилятор C++ можно попросить пренебречь достоверностью в ущерб скорости (fast math), но это скорее вольность языка.
Refridgerator
14.09.2023 11:23Это вы ещё си++ не смотрели, у которого в настройках компилятора аж целых 3 (три) режима работы с вещественными числами.
Vglk Автор
14.09.2023 11:23+1Смотрел, просто решил, что писать ещё и об этом - совсем статью раздувать. Там, помнится, можно указывать, использовать SSE или FPU, а ещё есть опция, которая позволяет компилятору ради скорости максимально упрощать выражения. С пометкой, что такое упрощение иногда может приводить к неправильным результатам. В частности, код для расчёта машинного эпсилон с этой опцией выдаёт 0, так как компилятор решает, что в выражении 1 + d/2 > 1 слева и справа можно сократить единицу, и упрощает его до d/2 > 0, а потом решает, что при сравнении с нулём на 2 делить не обязательно, и выражение превращается в d>0. В итоге деление идёт до тех пор, пока переменная не обратится в ноль.
adeshere
14.09.2023 11:23...есть опция, которая позволяет компилятору ради скорости максимально упрощать выражения (...) код для расчёта машинного эпсилон с этой опцией выдаёт 0
Это в не-интеловском компиляторе? Я проверял в интел-фортране в режиме Floating Point Model = Fast; Floating Point Speculation = Fast (вроде бы это аналог?) с оптимизацией на максимальную скорость. Там этот код корректно работает. Хотя проще наверно все-таки встроенную функцию вызвать:
r=EPSILON(r)
Конечно это фортран, но по идее, интеловский С++ должен так же себя вести?
Vglk Автор
14.09.2023 11:23+1Это я писал про компилятор из MS Visual Studio. Там есть опция /fp:fast. В справке про неё написано следующее:
Параметр
/fp:fast
позволяет компилятору изменять порядок, объединять или упрощать операции с плавающей запятой, чтобы оптимизировать код с плавающей запятой для скорости и пространства. Компилятор может пропускать округление при инструкциях присваивания, typecasts или вызовах функций. Он может переупорядочение операций или алгебраические преобразования, например, с помощью ассоциативных и распределительных законов. Он может изменить порядок кода, даже если такие преобразования приводят к заметно отличающимся поведением округления. Из-за этой расширенной оптимизации результат некоторых вычислений с плавающей запятой может отличаться от результатов, созданных другими/fp
вариантами. Специальные значения (NaN, +infinity, -infinity, -0,0) не могут распространяться или вести себя строго в соответствии со стандартом IEEE-754. Сокращения с плавающей запятой могут создаваться в ./fp:fast
Компилятор по-прежнему привязан к базовой архитектуре в/fp:fast
, и с помощью параметра могут быть доступны дополнительные оптимизации/arch
.В
/fp:fast
компилятор создает код, предназначенный для выполнения в среде с плавающей запятой по умолчанию, и предполагает, что среда с плавающей запятой не используется или не изменяется во время выполнения. То есть предполагается, что код: оставляет исключения с плавающей запятой маской, не считывает и не записывает регистры состояния с плавающей запятой и не изменяет режимы округления./fp:fast
предназначен для программ, которые не требуют строгого порядка исходного кода и округления выражений с плавающей запятой и не используют стандартные правила для обработки специальных значений, таких какNaN
. Если код с плавающей запятой требует сохранения порядка и округления исходного кода или использует стандартное поведение специальных значений, используйте/fp:precise
. Если код обращается к среде с плавающей запятой или изменяет ее для изменения режимов округления, распаковки исключений с плавающей запятой или проверка состояния с плавающей запятой, используйте ./fp:strict
Timofeuz
14.09.2023 11:23+1Значение типа decimal занимает в памяти 8 байт – 64 бита. Из них 1 бит отводится на знак, 96 – на целочисленное базовое значение B, 8 бит отводится на степень делителя F, 23 бита не используются.
Наверное всё-таки 16 байт -128 бит, или четыре int.
sergej0503
14.09.2023 11:23-1это что за платформа такая .Net? Программирую еще с начиная с .NET Framework 1.1 и не слышал о такой классификации, есть .NET Framework и .Net Core, а .Net это что за зверь или автору лень 5 дополнительных символов писать?
Vglk Автор
14.09.2023 11:23+3.NET Core уже давно официально переименован в просто .NET, ещё с 5-й версии. Вот пруф: https://learn.microsoft.com/ru-ru/lifecycle/products/microsoft-net-and-net-core
vivlav
Доброго дня, мне кажется есть небольшая путаница при разборе таких чисел в 10-м представлении таких чисел. (Разберу для double)
Рассмотрим ряд для 52 бит в двоичной системе:
0.1, 0.11, 0.111 ... 0.111..11{52 еденицы} им соответствуют десятичные:
0.5, 0.75, 0.875 ... 0,999..75{тоже 52 цифры!!!} то есть каждый бит в двоичной дроби - добавляет разряд в десятичной системе. И здесь и есть точное равенство двоичного представления.
Но представление в десятичное системе (как в статье) - обрежет и округлит их до 16 цифр в десятичном представлении. Это мы и видим во всем операциях-представлениям в статье все выводы в консоль - это десятичные представления округленные.
Точность в 16 знаков говорит о том, что изменяя любые цифры в числе выше (0,999..75) после 16-го знака в двоичной системе, нам будет давать то же самое число в переводе в двоичную систему. (На примере из статьи - можно сколько угодно добавлять знаков к десятичному представлению после 16-го символа). Поэтому для все конвертеры в десятичную строку и занимаются округлением.
Если отметить, как в статье, 0.1 в десятичной - бесконечная дробь в двоичной системе, то понятно почему 0.1f != 0.1d. (т.к. 0.00011001100..001{24 цифры для float} != 0.00011001100..001 {52 цифры для double}).
Может немного муторно написал, но хотел пояснить вот что:
0.1f - после перевода компилятором в двоичную будет выглядеть так:
0.000110011001100110011001101. Оно же в 10-й (в точности):
0.100000001490116119384765625
Представим, как это делает ToString():
0.1 - это округление для float
0.10000000149011612 - округление для double, хотя число одно и тоже в двоичном виде и целиком влезает во float и double. (можно проверить (double)(0.1f) == 0.1f).
Vglk Автор
Я, видимо, зря не заострил внимание на этом моменте в статье. Считать таким образом десятичные цифры бессмысленно, потому что имеют значение только те цифры, которые находятся в пределах погрешности, цифры, лежащие за погрешностью, значения не имеют. Когда, например, пишут, что float имеет точность 6-9 десятичных цифр, это значит, что от 6 до 9 цифр будут точно правильными, а дальше как повезёт. Как, например, для значения 0.1 отклонение начинается с 9-й значащей цифры. Для отдельных значений, являющихся суммами степеней двойки, точность может оказаться существенно выше и даже, как вы правильно заметили, дойти до 52 цифр. Но это - выбросы в статистике, для всех остальных чисел точность будет намного меньше.
vivlav
Понятно, что не имеет смысла (это только строго математически важно). Я скорее для понимания хотел добавить, что то, что выглядит разным в десятичном представлении, может быть одинаковым в двоичном.
double d = 0.1f;
float f = 0.1f;
Assert.True(d == f);//значения равны
Assert.False(d.ToString() == f.ToString());//представления не равны