Пользователь
0,0
рейтинг
1 февраля 2011 в 14:27

Разработка → Что нужно знать про арифметику с плавающей запятой из песочницы

C++*


В далекие времена, для IT-индустрии это 70-е годы прошлого века, ученые-математики (так раньше назывались программисты) сражались как Дон-Кихоты в неравном бою с компьютерами, которые тогда были размером с маленькие ветряные мельницы. Задачи ставились серьезные: поиск вражеских подлодок в океане по снимкам с орбиты, расчет баллистики ракет дальнего действия, и прочее. Для их решения компьютер должен оперировать действительными числами, которых, как известно, континуум, тогда как память конечна. Поэтому приходится отображать этот континуум на конечное множество нулей и единиц. В поисках компромисса между скоростью, размером и точностью представления ученые предложили числа с плавающей запятой (или плавающей точкой, если по-буржуйски).

Арифметика с плавающей запятой почему-то считается экзотической областью компьютерных наук, учитывая, что соответствующие типы данных присутствуют в каждом языке программирования. Я сам, если честно, никогда не придавал особого значения компьютерной арифметике, пока решая одну и ту же задачу на CPU и GPU получил разный результат. Оказалось, что в потайных углах этой области скрываются очень любопытные и странные явления: некоммутативность и неассоциативность арифметических операций, ноль со знаком, разность неравных чисел дает ноль, и прочее. Корни этого айсберга уходят глубоко в математику, а я под катом постараюсь обрисовать лишь то, что лежит на поверхности.

1. Основы


Множество целых чисел бесконечно, но мы всегда можем подобрать такое число бит, чтобы представить любое целое число, возникающее при решении конкретной задачи. Множество действительных чисел не только бесконечно, но еще и непрерывно, поэтому, сколько бы мы не взяли бит, мы неизбежно столкнемся с числами, которые не имеют точного представления. Числа с плавающей запятой — один из возможных способов предсталения действительных чисел, который является компромиссом между точностью и диапазоном принимаемых значений.

Число с плавающей запятой состоит из набора отдельных разрядов, условно разделенных на знак, экспоненту порядок и мантиссу. Порядок и мантисса — целые числа, которые вместе со знаком дают представление числа с плавающей запятой в следующем виде:



Математически это записывается так:

(-1)s × M × BE, где s — знак, B-основание, E — порядок, а M — мантисса.

Основание определяет систему счисления разрядов. Математически доказано, что числа с плавающей запятой с базой B=2 (двоичное представление) наиболее устойчивы к ошибкам округления, поэтому на практике встречаются только базы 2 и, реже, 10. Для дальнейшего изложения будем всегда полагать B=2, и формула числа с плавающей запятой будет иметь вид:

(-1)s × M × 2E

Что такое мантисса и порядок? Мантисса – это целое число фиксированной длины, которое представляет старшие разряды действительного числа. Допустим наша мантисса состоит из трех бит (|M|=3). Возьмем, например, число «5», которое в двоичной системе будет равно 1012. Старший бит соответствует 22=4, средний (который у нас равен нулю) 21=2, а младший 20=1. Порядок – это степень базы (двойки) старшего разряда. В нашем случае E=2. Такие числа удобно записывать в так называемом «научном» стандартном виде, например «1.01e+2». Сразу видно, что мантисса состоит из трех знаков, а порядок равен двум.

Допустим мы хотим получить дробное число, используя те же 3 бита мантиссы. Мы можем это сделать, если возьмем, скажем, E=1. Тогда наше число будет равно

1,01e+1 = 1×21+0×20+1×2-1=2+0,5=2,5

Здесь, поскольку E=1, степень двойки первого разряда (который идет перед запятой), равна «1». Два других разряда, расположенных правее (после запятой), обеспечивают вклад 2E-1 и 2E-2 (20 и 2-1 соответственно). Очевидно, что регулируя E одно и то же число можно представить по-разному. Рассмотрим пример с длиной мантиссы |M|=4. Число «2» можно представить в следующем виде:

2 = 10 (в двоичной системе) = 1.000e+1 = 0.100e+2 = 0.010e+3. (E=1, E=2, E=3 соответственно)

Обратите внимание, что одно и то же число имеет несколько представлений. Это не удобно для оборудования, т.к. нужно учитывать множественность представлния при сравнении чисел и при выполнении над ними арифметических операций. Кроме того, это не экономично, поскольку число представлений — конечное, а повторения уменьшают множество чисел, которые вообще могут быть представлены. Поэтому уже в самых первых машинах начали использовать трюк, делая первый бит мантиссы всегда положительным. Такое предаставление назвали нормализованным.



Это экономит один бит, так как неявную единицу не нужно хранить в памяти, и обеспечивает уникальность представления числа. В нашем примере «2» имеет единственное нормализованное представление («1.000e+1»), а мантисса хранится в памяти как «000», т.к. старшая единица подразумевается неявно. Но в нормализованном представлении чисел возникает новая проблема — в такой форме невозможно представить ноль.

Строго говоря, нормализованное число имеет следующий вид:

(-1)s × 1.M × 2E.

Качество решения задач во многом зависит от выбора представления чисел с плавающей запятой. Мы плавно подошли к проблеме стандартизации такого представления.

2. Немного истории


В 60-е и 70-е годы не было единого стандарта представления чисел с плавающей запятой, способов округления, арифметических операций. В результате программы были крайне не портабельны. Но еще большей проблемой было то, что у разных компьютеров были свои «странности» и их нужно было знать и учитывать в программе. Например, разница двух не равных чисел возвращала ноль. В результате выражения «X=Y» и «X-Y=0» вступали в противоречие. Умельцы обходили эту проблему очень хитрыми трюками, например, делали присваивание «X=(X-X)+X» перед операциями умножения и деления, чтобы избежать проблем.

Инициатива создать единый стандарт для представления чисел с плавающей запятой подозрительно совпала с попытками в 1976 году компанией Intel разработать «лучшую» арифметику для новых сопроцессоров к 8086 и i432. За разработку взялись ученые киты в этой области, проф. Джон Палмер и Уильям Кэхэн. Последний в своем интервью высказал мнение, что серьезность, с которой Intel разрабатывала свою арифметику, заставила другие компании объединиться и начать процесс стандартизации.

Все были настроены серьезно, ведь очень выгодно продвинуть свою архитектуру и сделать ее стандартной. Свои предложения представили компании DEC, National Superconductor, Zilog, Motorola. Производители мейнфреймов Cray и IBM наблюдали со стороны. Компания Intel, разумеется, тоже представила свою новую арифметику. Авторами предложенной спецификации стали Уильям Кэхэн, Джероми Кунен и Гарольд Стоун и их предложение сразу прозвали «K-C-S».

Практически сразу же были отброшены все предложения, кроме двух: VAX от DEC и «K-C-S» от Intel. Спецификация VAX была значительно проще, уже была реализована в компьютерах PDP-11, и было понятно, как на ней получить максимальную производительность. С другой стороны в «K-C-S» содержалось много полезной функциональности, такой как «специальные» и «денормализованные» числа (подробности ниже).

В «K-C-S» все арифметические алгоритмы заданы строго и требуется, чтобы в реализации результат с ними совпадал. Это позволяет выводить строгие выкладки в рамках этой спецификации. Если раньше математик решал задачу численными методами и доказывал свойства решения, не было никакой гарантии, что эти свойства сохранятся в программе. Строгость арифметики «K-C-S» сделала возможным доказательство теорем, опираясь на арифметику с плавающей запятой.

Компания DEC сделала все, чтобы ее спецификацию сделали стандартом. Она даже заручилась поддержкой некоторых авторитетных ученых в том, что арифметика «K-C-S» в принципе не может достигнуть такой же производительности, как у DEC. Ирония в том, что Intel знала, как сделать свою спецификацию такой же производительной, но эти хитрости были коммерческой тайной. Если бы Intel не уступила и не открыла часть секретов, она бы не смогла сдержать натиск DEC.

Подробнее о баталиях при стандартизации смотрите в интервью профессора Кэхэна, а мы рассмотрим, как выглядит представление чисел с плавающей запятой сейчас.

3. Представление чисел с плавающей запятой сегодня


Разработчики «K-C-S» победили и теперь их детище воплотилось в стандарт IEEE754. Числа с плавающей запятой в нем представлены в виде знака (s), мантиссы (M) и порядка (E) следующим образом:

(-1)s × 1.M × 2E

Замечание. В новом стандарте IEE754-2008 кроме чисел с основанием 2 присутствуют числа с основанием 10, так называемые десятичные (decimal) числа с плавающей запятой.

Чтобы не загромождать читателя чрезмерной информацией, которую можно найти в Википедии, рассмотрим только один тип данных, с одинарной точностью (float). Числа с половинной, двойной и расширенной точностью обладают теми же особенностями, но имеют другой диапазон порядка и мантиссы. В числах одинарной точности (float/single) порядок состоит из 8 бит, а мантисса – из 23. Эффективный порядок определяется как E-127. Например, число 0,15625 будет записано в памяти как


Рисунок взят из Википедии

В этом примере:
  • Знак s=0 (положительное число)
  • Порядок E=011111002-12710 = -3
  • Мантисса M = 1.012 (первая единица не явная)
  • В результате наше число F = 1.012e-3 = 2-3+2-5 = 0,125 + 0,03125 = 0,15625

Чуть более подробное объяснение
Здесь мы имеем дело с двоичным представлением числа «101» со сдвигом запятой на несколько разрядов влево. 1,01 — это двоичное представление, означающее 1×20 + 0×2-1 + 1×2-2. Сдвинув запятую на три позиции влево получим 1,01e-3 = 1×2-3 + 0×2-4 + 1×2-5 = 1×0,125 + 0×0,0625 + 1×0,03125 = 0,125 + 0,03125 = 0,15625.


3.1 Специальные числа: ноль, бесконечность и неопределенность

В IEEE754 число «0» представляется значением с порядком, равным E=Emin-1 (для single это -127) и нулевой мантиссой. Введение нуля как самостоятельного числа (т.к. в нормализованном представлении нельзя представить ноль) позволило избежать многих странностей в арифметике. И хоть операции с нулем нужно обрабатывать отдельно, обычно они выполняются быстрее, чем с обычными числами.

Также в IEEE754 предусмотрено представление для специальных чисел, работа с которыми вызывает исключение. К таким числам относится бесконечность (±∞) и неопределенность (NaN). Эти числа позволяет вернуть адекватное значение при переполнении. Бесконечности представлены как числа с порядком E=Emax+1 и нулевой мантиссой. Получить бесконечность можно при переполнении и при делении ненулевого числа на ноль. Бесконечность при делении разработчики определили исходя из существования пределов, когда делимое и делитель стремиться к какому-то числу. Соответственно, c/0==±∞ (например, 3/0=+∞, а -3/0=-∞), так как если делимое стремиться к константе, а делитель к нулю, предел равен бесконечности. При 0/0 предел не существует, поэтому результатом будет неопределенность.

Неопределенность или NaN (от not a number) – это представление, придуманное для того, чтобы арифметическая операция могла всегда вернуть какое-то не бессмысленное значение. В IEEE754 NaN представлен как число, в котором E=Emax+1, а мантисса не нулевая. Любая операция с NaN возвращает NaN. При желании в мантиссу можно записывать информацию, которую программа сможет интерпретировать. Стандартом это не оговорено и мантисса чаще всего игнорируется.

Как можно получить NaN? Одним из следующих способов:
  • ∞+(- ∞)
  • 0 × ∞
  • 0/0, ∞/∞
  • sqrt(x), где x<0

По определению NaN ≠ NaN, поэтому, для проверки значения переменной нужно просто сравнить ее с собой.

Зачем нулю знак (или +0 vs -0)

Любознательный читатель вероятно уже замелил заметил, что в описанном представлении чисел с плавающей запятой существует два нуля, которые отличаются только знаком. Так, 3·(+0)=+0, а 3·(-0)=-0. Но при сравнении +0=-0. В стандарте знак сохранили умышленно, чтобы выражения, которые в результате переполнения или потери значимости превращаются в бесконечность или в ноль, при умножении и делении все же могли представить максимально корректный результат. Например, если бы у нуля не было знака, выражение 1/(1/x)=x не выполнялось бы верно при x=±∞, так как 1/∞ и 1/-∞ равны 0.

Еще один пример:
(+∞/0) + ∞ = +∞, тогда как (+∞/-0) +∞ = NaN

Чем бесконечность в данном случае лучше, чем NaN? Тем, что если в арифметическом выражении появился NaN, результатом всего выражения всегда будет NaN. Если же в выражении встретилась бесконечность, то результатом может быть ноль, бесконечность или обычное число с плавающей запятой. Например, 1/∞=0.

3.3 Денормализованные числа

Что такое субнормальные денормализованные (subnormal) числа рассмотрим на простом примере. Пусть имеем нормализованное представление с длиной мантиссы |M|=2 бита (+ один бит нормализации) и диапазоном значений порядка -1≤E≤2. В этом случае получим 16 чисел:



Крупными штрихами показаны числа с мантиссой, равной 1,00. Видно, что расстояние от нуля до ближайшего числа (0 — 0,5) больше, чем от этого числа к следующему (0,5 — 0,625). Это значит, что разница двух любых чисел от 0,5 до 1 даст 0, даже если эти числа не равны. Что еще хуже, в пропасть между 0,5 и 0 попадает разница чисел, больших 1. Например, «1,5-1,25=0» (см. картинку).

В «околонулевую яму» подпадает не каждая программа. Согласно статистике 70-х годов в среднем каждый компьютер сталкивался с такой проблемой один раз в месяц. Учитывая, что компьютеры приобретали массовость, разработчики «K-C-S» посчитали эту проблему достаточно серьезной, чтобы решать ее на аппаратном уровне. Предложенное ими решение состояло в следующем. Мы знаем, что при E=Emin-1 (для float это «-127») и нулевой мантиссе число считается равным нулю. Если же мантисса не нулевая, то число считается не нулевым, его порядок полагается E=Emin, причем неявный старший бит мантиссы полагается равным нулю. Такие числа называются денормализованными.

Строго говодя, числа с плавающей запятой теперь имеют вид:

(-1)s × 1.M × 2E, если Emin≤E≤Emax (нормализованные числа)

(-1)s × 0.M × 2Emin, если E=Emin-1. (денормализованные числа)

Вернемся к примеру. Наш Emin=-1. Введем новое значение порядка, E=-2, при котором числа являются денормализованными. В результате получаем новое представление чисел:



Интервал от 0 до 0,5 заполняют денормализованные числа, что дает возможность не проваливаться в 0 рассмотренных выше примерах (0,5-0,25 и 1,5-1,25). Это сделало представление более устойчиво к ошибкам округления для чисел, близких к нулю.

Но роскошь использования денормализованного представления чисел в процессоре не дается бесплатно. Из-за того, что такие числа нужно обрабатывать по-другому во всех арифметических операциях, трудно сделать работу в такой арифметике эффективной. Это накладывает дополнительные сложности при реализации АЛУ в процессоре. И хоть денормализованные числа очень полезны, они не являются панацеей и за округлением до нуля все равно нужно следить. Поэтому эта функциональность стала камнем преткновения при разработке стандарта и встретила самое сильное сопротивление.

3.4 Очередность чисел в IEEE754

Одна из удивительных особенностей представления чисел в формате IEEE754 состоит в том, что порядок и мантисса расположены друг за другом таким образом, что вместе образуют последовательность целых чисел {n} для которых выполняется:

n<n+1 ⇒ F(n) < F(n+1), где F(n) – число с плавающей запятой, образованное от целого n, разбиением его битов на порядок и мантиссу.

Поэтому если взять положительное число с плавающей запятой, преобразовать его к целому, прибавить «1», мы получим следующее число, которое представимо в этой арифметике. На Си это можно сделать так:

float a=0.5;
int n = *((int*) &a);
float b = *((float*) &(++n));
printf("После %e следующее число: %e, разница (%e)\n", a, b, b-a);

Этот код будет работать только на архитектуре с 32-битным int.

4. Подводные камни в арифметике с плавающей запятой


Теперь – к практике. Рассмотрим особенности арифметики с плавающей запятой, к которым нужно проявить особую осторожность при программировании.

4.1 Округление

С ошибками из-за погрешностей округления в современной арифметике с плавающей запятой встретиться сложно, особенно если использовать двойную точность. Правило округления в стандарте IEEE754 говорит о том, что результат любой арифметической операции должен быть таким, как если бы он был выполнен над точными значениями и округлен до ближайшего числа, представимого в этом формате. Это требует от АЛУ дополнительных усилий и некоторые опции компилятора (такие как «-ffast-math» в gcc) могут отключить такое поведение. Особенности округления в IEEE754:
  • Округление до ближайшего в стандарте сделано не так как мы привыкли. Математически показано, что если 0,5 округлять до 1 (в большую сторону), то существует набор операций, при которых ошибка округления будет возрастать до бесконечности. Поэтому в IEEE754 применяется правило округления до четного. Так, 12,5 будет округлено до 12, а 13,5 – до 14.
  • Самая опасная операция с точки зрения округления в арифметике с плавающей запятой — это вычитание. При вычитании близких чисел значимые разряды могут потеряться, что
    может в разы увеличить относительную погрешность.
  • Для многих широко распространенных математических формул математики разработали специальную форму, которая позволяет значительно уменьшить погрешность при округлении. Например, расчет формулы «x2-y2» лучше вычислять используя формулу «(x-y)(x+y)».

4.2 Неассоциативность арифметических операций

В арифметике с плавающей запятой правило (a*b)*c = a*(b*c) не выполняется для любых арифметических операций. Например,

(1020+1)-1020=0 ≠ (1020-1020)+1=1

Допустим у нас есть программа суммирования чисел.

double s = 0.0;
for (int i=0; i<n; i++) s = s + t[i];

Некоторые компиляторы по умолчанию могут переписать код для использования нескольких АЛУ одновременно (будем считать, что n делится на 2):

double sa[2], s; 
sa[0]=sa[1]=0.0;
for (int i=0; i<n/2; i++) {
    sa[0]=sa[0]+t[i*2+0];
    sa[1]=sa[1]+t[i*2+1];
}
S=sa[0]+sa[1];

Так как операции суммирования не ассоциативны, эти две программы могут выдать различный результат.

4.3 Числовые константы

Помните, что не все десятичные числа имеют двоичное представление с плавающей запятой. Например, число «0,2» будет представлено как «0,200000003» в одинарной точности. Соответственно, «0,2 + 0,2 ≈ 0,4». Абсолютная погрешность в отдельном
случае может и не высока, но если использовать такую константу в цикле, можем получить накопленную погрешность.

4.4 Выбор минимума из двух значений

Допустим из двух значений нам нужно выбрать минимальное. В Си это можно сделать одним из следующих способов:
  1. x < y? x: y
  2. x <= y? x: y
  3. x > y? y: x
  4. x >= y? y: x

Часто компилятор считает их эквивалентными и всегда использует первый вариант, так как он выполняется за одну инструкцию процессора. Но если мы учтем ±0 и NaN, эти операции никак не эквивалентны:
x y x < y? x: y x <= y? x: y x > y? y: x x >= y? y: x
+0 -0 -0 +0 +0 -0
NaN 1 1 1 NaN NaN

4.5 Сравнение чисел

Очень распространенная ошибка при работе с float-ами возникает при проверке на равенство. Например,

float fValue = 0.2;
if (fValue == 0.2) DoStuff();

Ошибка здесь, во-первых, в том, что 0,2 не имеет точного двоичного представления, а во-вторых 0,2 – это константа двойной точности, а переменная fValue – одинарной, и никакой гарантии о поведении этого сравнения нет.

Лучший, но все равно ошибочный способ, это сравнивать разницу с допустимой абсолютной погрешностью:

if (fabs(fValue – fExpected) < 0.0001) DoStuff(); // fValue=fExpected?


Недостаток такого подхода в том, что погрешность представления числа увеличивается с ростом самого этого числа. Так, если программа ожидает «10000», то приведенное равенство не будет выполняться для ближайшего соседнего числа (10000,000977). Это особенно актуально, если в программе имеется преобразование из одинарной точности в двойную.

Выбрать правильную процедуру сравнения сложно и заинтересованных читателей я отсылаю к статье Брюса Доусона. В ней предлагается сравнивать числа с плавающей запятой преобразованием к целочисленной переменной. Это — лучший, хотя и не портабельный способ:

bool AlmostEqual2sComplement(float A, float B, int maxUlps) {
    // maxUlps не должен быть отрицательным и не слишком большим, чтобы
    // NaN не был равен ни одному числу
    assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
    int aInt = *(int*)&A;
    // Уберем знак в aInt, если есть, чтобы получить правильно упорядоченную последовательность
    if (aInt < 0) aInt = 0x80000000 - aInt;
    //aInt &= 0x7fffffff; //(см. комментарий пользователя Vayun)
    // Аналогично для bInt
    int bInt = *(int*)&B;
    if (bInt < 0) bInt = 0x80000000 - bInt;
    /*aInt &= 0x7fffffff;*/
    unsigned int intDiff = abs(aInt - bInt); /*(см. комментарий пользователя Vayun)*/
    if (intDiff <= maxUlps)
        return true;
    return false;
}


В этой программе maxUlps (от Units-In-Last-Place) – это максимальное количество чисел с плавающей запятой, которое может лежать между проверяемым и ожидаемым значением. Другой смысл этой переменной – это количество двоичных разрядов (начиная с младшего) в сравниваемых числах разрешается упустить. Например, maxUlps=16, означает, что младшие 4 бита (log216) могут не совпадать, а числа все равно будут считаться равными. При этом, при сравнении с числом 10000 абсолютная погрешность будет равна 0,0146, а при сравнении с 0.001, погрешность будет менее 0.00000001 (10-8).

5. Проверка полноты поддержки IEE754


Думаете, что если процессоры полностью соответствуют стандарту IEEE754, то любая программа, использующая стандартные типы данных (такие как float/double в Си), будет выдавать один и тот же результат на разных компьютерах? Ошибаетесь. На портабельность и соответствие стандарту влияет компилятор и опции оптимизации. Уильям Кэхэн написал программу на Си (есть версия и для Фортрана), которая позволяет проверить удовлетворяет ли связка «архитектура+компилятор+опции» IEEE754. Называется она «Floating point paranoia» и ее исходные тексты доступны для скачивания. Аналогичная программа доступна для GPU. Так, например, компилятор Intel (icc) по умолчанию использует «расслабленную» модель IEEE754, и в результате не все тесты выполняются. Опция «-fp-model precise» позволяет компилировать программу с точным соответствием стандарту. В компиляторе GCC есть опция «-ffast-math», использование которой приводит к несоответствию IEEE754.

Заключение


Напоследок поучительная история. Когда я работал над тестовым проектом на GPU, у меня была последовательная и параллельная версия одной программы. Сравнив время выполнения, я был очень обрадован, так как получил ускорение в 300 раз. Но позже оказалось, что вычисления на GPU «разваливались» и обращались в NaN, а работа с ними в GPU была быстрее, чем с обычными числами. Интересно было другое — одна и та же программа на эмуляторе GPU (на CPU) выдавала корректный результат, а на самом GPU – нет. Позже оказалось, что проблема была в том, что этот GPU не поддерживал полностью стандарт IEEE754 и прямой подход не сработал.

Сейчас арифметика с плавающей запятой почти совершенна. Практически всегда наивный подход сработает, и программа, не учитывающая все ее особенности, выдаст правильный результат, а описанные подводные камни касаются только экзотических случаев. Но нужно всегда оставаться бдительным: в таком вопросе как компьютерная математика легко наступить на грабли.

P.S. Спасибо пользователю uqlock за важное замечание. Деньги нельзя хранить в виде числа с плавающей запятой, т.к. в этом случае нельзя выделить значимые разряды. Если в языке программирования нет типов данных с фиксированной запятой, можно выйти из положения и хранить деньги в виде целого числа, подразумевая копейки (иногда доли копеек).

P.P.S. Благодарность за найденные опечатки и ошибки пользователям: gribozavr, kurokikaze, Cenness, TheShock, perl_demon, GordTremor, fader44, DraculaDis, icc, f0rbidik, Harkonnen, AlexanderYastrebov, Vayun, EvilsInterrupt!

Литература

  1. Интервью с Уильамом Кэхэном про становление стандарта IEE754.
  2. What Every Computer Scientist Should Know About Floating-Point Arithmetic, David Goldberg — книга с математическими выкладками.
  3. Comparing floating point numbers, Bruce Dowson.
Руслан Ющенко @yruslan
карма
106,7
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама

Самое читаемое Разработка

Комментарии (74)

  • +7
    В принципе, достаточно доходчиво объяснено. А то оооочень много личностей даже понятия не имеют что там внутри double. Буду теперь новичкам этот пост показывать. Спасибо.
    • 0
      Очень интересный пост, очень понятно объясняется и главное по полочкам.

      Но все равно силился час пока прочел и смог запомнить хотя бы часть.

      Огромное спасибо.
      • +1
        Огромное спасибо за пост! Пару поправок:
        1. портабельность
        2. int должен быть не только 32-битным, но и little-endian (x86, к примеру)
        • 0
          Про порядок байт — под вопросом
          • 0
            Смотря оговорена ли endianness для float'ов в IEEE или же наследует int'овую
        • 0
          Портабельность → переносимость
  • +1
    Мне всегда казалось, что Floating Point всеже подразумевает плавающую точку.
    • 0
      В смысле термин буржуйский и перевод нужен нормальный, а то както в шок вводит немного.
      • +4
        Я консультировался с разработчиками стандартов в Украине. Оба перевода правильные, а в стандартах пишут через слеш (плавающая точка/запятая). Последний мне показался самобытнее и ближе к традиции.
        • –16
          лучше бы проконсультировался с правилами русского языка…

          • –5
            убогие каклы :-)
        • –1
          > Любознательный читатель вероятно уже замеЛил
        • 0
          >Последний мне показался самобытнее и ближе к традиции.

          Как по мне, так лучше бы «ближе к общепринятой терминологии» (а то «ближе к традиции» порой приводит к «блогозаписи» и прочим «междумордиям»).

          Но на то вы и автор, чтобы решать самостоятельно.
      • 0
        Ну это как бы не для лингвистов перевод. А для программистов, которые ещё со времён совка называем это арифметикой с плавающей запятой.
        Давайте ещё команды ассемблера на русский язык, корректно переведём.
        • НЛО прилетело и опубликовало эту надпись здесь
          • НЛО прилетело и опубликовало эту надпись здесь
      • НЛО прилетело и опубликовало эту надпись здесь
      • 0
        Нормальный перевод. У них запятая иногда — разделитель троек десятичных разрядов. У нас же запятая отделяет дробную часть. Так что 1,001 чревато мисспеллингом =)
        • 0
          Насколько мне известно, до 10 000 (10,000) числа пишутся слитно.
  • +4
    Статья — торт.
    • +6
      Надо на Хабре добавить мегатонны смайликов, как в некоторых развлекательных блогах, и вообще можно будет не думать, что написать в комментарии.
  • +7
    >компьютер должен оперировать рациональными числами, которых, как известно, континуум

    Извиняюсь за занудство, но рациональных чисел все-таки счетное число. Это поменьше чем континуум.
    • 0
      Спасибо! Уже исправил. Множество действительных чисел, конечно.
      • +1
        > Множество действительных чисел, конечно.

        [каламбур]Множество действительных чисел бесконечно[/каламбур]
        Велик и могуч русский язык :) Но как раз слово «конечно» в обиходе является синонимом «обозримо мозгом, имеет начало и конец» и поэтому будет антонимом для «несчетно, неисчислимо, бесконечно».
  • +21
    Самое главное не сказано в статье — нельзя использовать float/double для рассчетов с деньгами. Только числа с фиксированой запятой. Иначе на огруглениях можно много потерять и потом баланс не сойдется.
    • –1
      Согласен. Тем более, если мы будем считать миллиарды, можем потерять не только копейки, но и рубли.
      • 0
        на этом в Англии один подросток банк на несколько тысяч фунтов обставил
      • +1
        вообще-то double (64-bit) хранит 53 двоичных цифры (+ 11 цифр экспоненты + 1 цифру знака), а это чуть больше 19 десятичных цифр. Соответственно, при счёте миллиардов мы не только не потеряем ни копейки, но даже учтём миллиардные доли копейки
        • 0
          чёрт, соврал, речь про extended — 80-bit — именно такие регистры в математических сопроцессорах, наследниках i8087. Они помнят по 19 знаков.

          double, конечно, помнит чуть больше 15 знаков. Ну то есть, если считать миллиарды, то учтём только лишь стотысячные доли копейки, какая жалость
          • 0
            Все зависит от того, что мы с этими числами будем делать. Кроме округления при представлении источником ошибок являются арифметические операции, которые тоже вносят погрешность. Эффект получается кумулятивный.
          • +1
            >Ну то есть, если считать миллиарды, то учтём только лишь стотысячные доли копейки, какая жалость

            Что случится после 100 миллионов операций с погрешностью в стотысячные доли копеек каждая?
            • 0
              Во-первых, это зависит от рода операций, а во-вторых — погрешности не растут линейно.

              Если вы имеете ввиду сложение близких по порядку чисел, т.е. определённых с одинаковой абсолютной погрешностью, то абсолютная погрешность их суммы будет в sqrt(100 миллионов) раз больше.
              • 0
                >… во-вторых — погрешности не растут линейно.

                Конечно, не линейно. Иногда они имеют свойство перемножаться, т.е. расти экспоненциально :) Грубо говоря, не «100М * 0.0000001», а «1.0000001^100М», а это уже са-а-аффсем другой разговор получится.

                Вы, безусловно, правы, и всё зависит от рода операций, но я думаю, финансисты будут очень возражать, если им сказать «наша система в целом работает нормально, но на некоторых операциях может накапливаться погрешность» %)
                • 0
                  Хорошо, предположим что ваш тезис о недопустимости финансовых расчетов с помощью арифметики с плавающей точкой верен. А какие в таком случае альтернативы? Если имеется счет, на который была внесена сумма в ххх.xx рублей с процентной ставкой yy.yy% в месяц и накопительным процентом (я забыл как называется этот термин правильно) — то какая сумма должна быть на счете через z1 дней, если через z2 дней клиент захочет снять треть вклада? Посчитайте это без использования плавающей точки, перемножьте на количество клиентов в банке и сравните получившиеся погрешности. Так или иначе придется «вручную» совершать те же операции нормализации и округления до приемлемых порядков.
                  • +1
                    На самом деле главный затык даже не в округлении как таковом, а именно в плавающей точке и в ограниченности длины мантиссы. При добавлении малых величин к большим, младшие разряды могут банально потеряться, т.к. просто не влезут в мантиссу.

                    Простой пример: просуммировать числа вида 1/n при n от 1 до какого-нибудь K. При малых K нет никакой разницы, суммировать от меньшего к большему или от большего к меньшему. Но как только K cтановится достаточно большим, разница становится заметной на глаз.

                    Простейший тест:

                    #include <stdio.h>
                    
                    void main(int argc, char* argv[]) {
                    	float sum1 = 0.f;
                    	float sum2 = 0.f;
                    	int n = 20000;
                    	int i;
                    
                    	printf("%i numbers\n", n);
                    
                    	for ( i = 1; i <= n; ++i ) {
                    		sum1 += 1.f / i;
                    		sum2 += 1.f / (n + 1 - i);
                    	}
                    	printf("%f\n", sum1 - sum2);
                    }


                    даст разницу 0.000007 для 1К итераций, 0.000009 для 10К итераций и 0.000023 для 20К итераций. А казалось бы, «от перестановки мест слагаемых сумма не меняется»… :)

                    Происходит это вот почему:
                    При подсчёте sum1 в неё сразу кладётся 1, фиксируя экспоненту, после чего начинают добавляться всё более мелкие числа. В итоге последние из них банально «обрубаются» (если не полностью, то частично), т.к. попросту не влазят в мантиссу.
                    При подсчёте sum2, напротив, сначала берутся более мелкие числа, которые накапливаются постепенно, отбрасывая младшие разряды по одному по мере роста экспоненты.

                    Конечно, double немного «отодвигает» эту проблему, но в некоторых особо точных расчётах и его может не хватить (особенно с учётом накладывающейся погрешности округления).
                    • 0
                      > но в некоторых особо точных расчётах и его может не хватить

                      вот это я и пытаюсь спросить — каким образом вы собираетесь делать «особо точные расчеты», если не использовать те же алгоритмы с увеличенной разрядностью? ладно, умножение в столбик я представляю как реализовать. а деление и прочие синусы?
                      • 0
                        >вот это я и пытаюсь спросить — каким образом вы собираетесь делать «особо точные расчеты», если не использовать те же алгоритмы с увеличенной разрядностью?

                        Считать и целую, и дробную части как отдельные целые (сбрасывая определённый старший разряд из копеек в рубли переполнении). Тем самым можно и рубли посчитать до милионов, и копейки до миллионных долей, «мантисса» получается достаточно широкой, а точка зафиксирована (в простейшем случае — ровно на середине).
                    • 0
                      Бухгалтерские рассчёты традиционно организованы как операции с фиксированной запятой. Именно так:
                      При подсчёте sum1 в неё сразу кладётся 1, фиксируя экспоненту, после чего начинают добавляться всё более мелкие числа. В итоге последние из них банально «обрубаются» (если не полностью, то частично), т.к. попросту не влазят в мантиссу.

                      только лишь с той разницей, что достаточно мелкие числа вообще непредставимы при таком подходе. Фиксируя запятую, мы фиксируем экспоненту. Фиксированная запятая означает фиксированную абсолютную погрешность. Т.е. считаем мы всегда с точностью «до копейки», неважно, что именно считаем — копейки или миллиарды рублей.

                      Мы не приобретаем дополнительной точности при округлении. Напротив, для достаточно мелких (но ещё представимых) чисел округление приводит к катастрофическому росту относительной погрешности. Как вы понимаете, этого не происходит для чисел с плавающей запятой, кроме денормализованных, которые фактически есть числа с фиксированной запятой. Для них округление всегда происходит в далёком младшем разряде, относительная погрешность округления всегда не более 2^-52 (для double).

                      Для организации фиксированной запятой используется обычная целочисленная арифметика, которая гораздо проще, а значит надёжнее и быстрее. (Запятая фактически появляется только при выводе чисел на экран.)
                      Этого вполне достаточно для бухгалтерии. Там очень сильно ограничен диапазон чисел, с которыми они работают, и большой диапазон чисел, который даёт плавающая запятая, просто не нужен. А (не нужной им) точностью они готовы пожертвовать в пользу надёжности и скорости (вспомните, всё это образовалось полвека назад, когда компьютеры были большими).
    • +6
      Надо показать эту статью Навальному =)
    • +2
      Добавлю, что в .NET можно использовать decimal для рассчетов с деньгами
  • +2
    надо было еще про денормализованные числа написать
    • +1
      См. п.3.3. Вы правы, денормализованные (вместо субнормальные) числа — правильный термин. Исправил в статье. Спасибо!
  • 0
    Всегда было гораздо интереснее — как именно нужно изменить знакомые формулы, чтобы компьютер с плавающей запятой считал по ним точнее (ну хотя бы вообще получал хоть что-то похожее).

    А то есть вот такая формула для синуса (ряд Маклорена): x-x^3/3!+x^5/5!-…
    Хорошая формула. Радиус сходимости — бесконечность. Только вот если считать с её помощью на компьютере (скажем, заканчивать, когда очередное слагаемое меньше затребованной погрешности), чуть отойдёшь от нуля — и всё, получаются совершенно несуразные результаты…
    • 0
      Простого ответа нет — есть целая теория как делать формулы более устойчивыми к арифметике с плавающей запятой, основанная на анализе погрешности и алгебре. Несколько примеров с доказательствами есть в этой книжке.
      • 0
        Я понимаю, что «не так-то просто». Я как бы намекаю, как можно развивать тему, начатую этой статьей :)
    • 0
      Это наука называется «вычислительные методы» :-)
  • 0
    Уильям Кэхэн написал программу на Си (есть версия и для Фортрана), которая позволяет проверить удовлетворяет ли связка «архитектура+компилятор+опции» IEEE754.

    От людей, занимающихся мат моделированием на физ-техе, слышал, что Fortran используется для этих целей как раз потому, что на всех платформах выдаёт один и тот же результат для операций с плавающей точкой. Это его главное преимущество в глазах учёных перед другими языками. Однако из ваших слов следует, что это не так. Объясните, пожалуйста: меня обманули? Или чего-то недосказали?
    • +2
      быть не может. Фортран — компилируемый, значит в итоге получим машинный код, который выполняется на конкретной машине, у которой конкретная архитектура АЛУ, либо если оное отсутствует, какая-то конкретная библиотека эмуляции арфиметики. Вот и получили запвисимость от реализации.
    • +1
      Фортран используется исключительно по традиции:

      в университетах живёт, здравствует и преподаёт куча именитых аксакалов, которые слышать не хотят про C, Python и т. п. и библиотеки вычислений с произвольной точностью

      есть куча прекрасной учебной и справочной литературы по численным методам с примерами и даже готовыми текстами программ на Фортране. В качестве примера вспоминается двухтомник Хайрера, Нёрсета и Виннера про численное решение систем дифференциальных уравнений с упоминанием исходных текстов знаменитых и нашумевших в своё время программ Дормана-Принса

      есть куча самих этих программ и библиотек на Форнтране. Хотя некоторые библиотеки в принципе можно скомпоновать с модулями, написанными на других языках, есть и такие, которые не скомпонуешь
  • +1
    Простите за занудство, но у вас в табличке «x > y? y: x» повторяется дважды.
    Уж больно статья хорошая, хочется, чтобы даже досадных опечаток в ней не было. :)
    • 0
      Спасибо! Исправил.
  • 0
    Про округление: помимо описанного способа округления, используемого по умолчанию, стандарт IEEE 754-2008 предусматривает еще 4 режима округления: в сторону 0, +inf, -inf и до ближайшего в сторону 0.

    Также стандарт определяет два типа NaN: Quiet и Signaling. Первый молча передает значение NaN в последующие вычисления, в то время как операции над вторым приводят к исключению.

    Про денормализованные числа полезно знать, что работа с ними, даже если она реализована в железе, как правило, выполняется гораздо медленнее (в десятки раз), чем с нормализованными числами. Поэтому в приложениях, нечувствительных к точности, но где важно быстродействие (3D графика, обработка сигналов), зачастую используется режим flush-to-zero. Многие FPU такой режим поддерживают аппаратно.

    Это так, небольшие дополнения. Статья суперская!
  • +1
    Что такое субнормальные денормализованные (subnormal) числа рассмотрим на простом примере. Пусть имеем нормализованное представление с длиной мантиссы P=3 (бита) и диапазоном значений экспоненты -1≤E≤2. В этом случае получим 16 чисел:

    В вашем примере длина мантиссы 2 бита.

    float a=0.5;
    int n = *((int*) &a);
    float b = *((int*) &(++n));
    ...
    

    s/float b = *((int*) &(++n));/float b = *((float*) &(++n));/

        ...
        if (aInt < 0)
            aInt = 0x80000000 - aInt;
        ...
        if (bInt < 0)
            bInt = 0x80000000 - bInt;
        ...
    

    Для знаковых целых:

    0x80000000 <= aInt
    0x80000000 <= bInt
    

    Поэтому должна быть либо сумма:

    aInt = 0x80000000 + aInt;
    ...
    bInt = 0x80000000 + bInt;
    

    либо разность с обратным знаком:

    aInt = aInt - 0x80000000;
    ...
    bInt = bInt - 0x80000000;
    

    Ну, или:

        ...
        /*
        if (aInt < 0)
            aInt = 0x80000000 - aInt;
        ...
        */
        aInt &= 0x7fffffff;
        ...
        /*
        if (bInt < 0)
            bInt = 0x80000000 - bInt;
        ...
        */
        bInt &= 0x7fffffff;
    
    • 0
      Да, вы правы. Большое спасибо за замечания!
      • 0
        Вам спасибо, статья отличная!
      • +1
        статья полезная, но не стоило менять правильное на неправильное))

        В приведенной ссылке www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm, функция AlmostEqual2sComplement написана корректно.

        Там именно должно быть
        if (aInt < 0)
        aInt = 0x80000000 - aInt;


        Потому что порядок float совпадает с порядком sign-and-magnitude integer. В то время как процессоры оперируют two's complement integer. Чтобы найти расстояние между двумя float в ulp, нужно одно перевести в другое и это делается через 0x80000000 - aInt для отрицательных float.

        Обрезка знакового бита через &= 0x7fffffff; приведет к тому, что функция будет считать равными числа отличающиеся знаком (например +1 будет равен -1)

        Плюс, есть еще замечательные грабли, которые имеются и у Bruce Dawson.
        int intDiff должен быть unsigned int, тогда будет использоваться беззнаковое сравнение (по правилам сравнения С), иначе если в abs случится переполнение, её результат будет интерпретирован как отрицательный и if (intDiff <= maxUlps) вернет истину для совершенно разных чисел.
        • 0
          Большое спасибо за содержательную критику!

          Я внес исправления в программу, но теперь возникает вопрос с исключением. Не возникнет ли целочисленное переполнение в этом случае? И не приведет ли это к замедлению программы?
          Будет время, сосредоточусь и проверю. )
          • +1
            Целочисленные операции не вызывают исключения при переполнении (по крайней мере на архитектурах, где описываемый трюк работает) и соответственно не будет никакого замедления, просто результат обрежется.

            В нашем случае все еще проще, результат который возвращает abs может не поместиться в диапазон положительных значений int, но unsigned int точно хватит.
            (тк сумма или разность может быть в общем случае длиннее на 1 бит, например 0b1111+0b1111=0b11110)

            Вообще в x86(_64) операции сложения и вычитания одинаковы для знаковых и беззнаковых операндов (собственно поэтому там используется two's complement форма для отрицательных чисел — оптимизация).

            Но вот операции сравнения уже различны. Для одного и того же бинарного значения в регистре разные операции сравнения интерпретируют его как знаковое или беззнаковое. Какую операцию использовать мы указываем компилятору через тип переменной (кстати интересно почему abs возвращает int а не unsigned).

            ЗЫ все это достаточно сложно прочитать с листа. Легче всего написать тестовую программу. Для float можно сделать почти полный перебор всех возможных входных значений)_
  • 0
    А вот мы сейчас это все будем читать, долго и вдумчиво. И наверное даже отправим на Киндл, да. А автору поставим плюс в карму фигурально и буквально, ибо верной дорогой идет.
  • 0
    Отличный пост.
  • 0
    Сейчас при приеме вычислительных техзаданий появилась новая мода как Bitwise Reproducibility, это полная ж для оптимизирующих компиляторов, которые любят менять местами порядок выполнения операций.

    Кстати есть отдельный вопрос как ловить NaN'ы (мы их так и читаем нанЫ :)
    -Слышь, у меня нанЫ в программе повалили, шоделать?).
    В фортране есть ключи компилятора, которые позволяют вываливаться в эксепшн при возникновении NaN/Inf. Для gfort -ffpe-trap=invalid,zero,overflow, для ifort -fpe0, для C/C++ несколько сложнее, нужно менять код (в районе main()) и вызывать такую вот функцию:
    feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);

    Еще бы на целочисленную арифметику (переполнение) уметь генерировать эксепшны… Но это я уже мечтаю, да :)

    • 0
      зачем эксепшены? в x86 для этого carry-флаг есть.
        add    RAX,RBX
        jc     >label
      
      Надо понимать, что «сложение» и «сложение xx-битных чисел в компьютере» вообще говоря разные операции. Эксепшен — это «все, финиш, приехали». А переполнение при том же сложении — это заранее оговоренная в документации особенность работы конкретного процессора.
      • 0
        Окей, как заставить один из общеупотребимых компиляторов проверять этот флаг?
        А почему эксепшн — так если у меня индекс переполняется, то мне как раз и нужно «все, приехали» — выпадение в кор, с генерацией дампа, бэктрейса, чтобы подставив отладчик я за два прогона нашел в каком месте оно валится и какой локал у меня при этом.
    • 0
      NaN с любым логическим выражением вернет ложь. Отсюда,
      if (a!=a) printf("Оппа, NaN..."); 
      • 0
        не во всех языках есть NaN

        например перловый undef не есть NaN в чистом виде, там undef вполне определенное значение… (ой, а может оно тоже последние биты теряет… в ядро не лазил, сорри)

        а для полезности перловиков, да и не только их — используйте Big Int для денег и их копеек — оно разряды не теряет ибо их нет.
      • +2
        Небольшое уточнение. На самом деле операция с NaN не возвращает ложь. Это компилятор так трактует. Не могу сказать по поводу SSE, но на FPU операция сравнения возвращает 3 флага:
        www.felixcloutier.com/x86/FCOM:FCOMP:FCOMPP.html
        То есть когда мы пишем:
        if (a==a) printf(«не NaN»)
        мы можем получить 3 разных состояния в регистрах:
        C3 = 0; C2 = 0; C0 = 0 — значит числа не равны
        С3 = 1; C2 = 0; C0 = 0 — значит числа равны
        С3 = 1; C2 = 1; C0 = 1 — как минимум один из операндов NaN
        Т.е. это 3 разных результата, и обычно компиляторы третье состояние трактуют как False (несмотря на то что C3 = 1, что соответствует True для операции сравнения).
  • +1
    Я понимаю, топик про float, но, взглянув на картинку, не могу удержаться.

    CL-USER(1): (+ 1/5 1/5)
    2/5

    :P
  • 0
    Однозначно торт.
  • 0
    Какие преимущества дает смещение порядка?
    • 0
      Один дополнительный бит для поля степени. Иначе надо было бы выделять бит для знака степени.
  • 0
    Извините за тупость и глупый вопрос, но разложите кто-нибудь более детально формулу получения из 1.01e-3 числа 0,15625. Очень хочется разобраться во всем, но своих мозгов что-то не хватает…
    • +2
      Здесь мы имеем дело с двоичным представлением числа «101» со сдвигом запятой на несколько разрядов влево. 1,01 — это двоичное представление, означающее 1*20 + 0*2-1 + 1*2-2 (где * — умножение). Сдвинув запятую на три позиции влево получим «1,01e-3» = 1*2-3 + 0*2-4 + 1*2-5 = 1*0,125 + 0*0,0625 + 1*0,03125 = 0,125 + 0,03125 = 0,15625.
      • 0
        Огромное Вам спасибо! Теперь все прояснилось!

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.