vladec Oct 25 2011 at 23:12

Вычисления с фиксированной точкой. Основные принципы (ч.1)

11 min

90K

Введение или зачем этот топик

Читая Хабрахабр, я натолкнулся на два топика, «выводящие на чистую воду» вычисления с плавающей запятой.
В одном из них достаточно подробно и качественно дана выжимка из стандарта IEEE754 и основные проблемы при вычислениях с плавающей запятой, другой — короткий топик-заметка про то, что не все так хорошо при вычислениях на ПК. При этом даются рекомендации в случае, когда важна математическая точность результата, использовать целочисленные вычисления, «фиксировать запятую» или как минимум проверять результаты, выдаваемые платформой (компилятор + процессор).
Несмотря на то, что советы дельные, понять, как использовать целочисленные вычисления там, где до этого была плавающая запятая, не так просто, особенно без математической подготовки. Достаточно занимательна в этом смысле попытка одного из «хабровчан» разобраться с фиксированной точкой методом экспериментов.
Данный топик — краткое введение, которое должно дать представление о вычислениях с фиксированной точкой. Математика в данной статье не должна никого напугать — все очень примитивно. Сразу прошу простить: среди моих знакомых устоявшимся выражением является именно «фиксированная точка» (от англ., fixed-point), а не «запятая», поэтому я буду придерживаться именно этого термина.

Еще раз о мантиссе и экспоненте

В вычислительной математике дробные значения представляют в виде пары целых чисел (n, e): мантиссы и экспоненты (по-русски более верно «показателя степени», но для краткости и по привычке буду в дальнейшем употреблять именно слово «экспонента»). Пара представляет дробное число в виде n * 2^-e.
Экспоненту можно рассматривать как количество цифр перед запятой, отделяющей дробную часть числа.
Если экспонента переменная, записываемая в регистр и неизвестная при компиляции, (n, e) называют числом с плавающей запятой. Если экспонента известна заранее, (n, e) называют числом с фиксированной точкой. Числа с фиксированной точкой могут записываться в обыкновенные целочисленные переменные (регистры) путем сохранения только мантиссы. Экспоненту обычно обозначают буквой q. Так что, встретив в комментарии к переменной что-нибудь в духе "q15 multiplier", следует рассматривать эту переменную как число с фиксированной точкой и экспонентой, равной 15. Впрочем, я еще вернусь к вопросу нотации, встречающейся в различных исходниках и статьях.

Вычисления

Итак, мы разобрались с тем, что при работе с фиксированной точкой экспонента нигде не записывается и держится «в уме».
Как же производить расчеты? Вычислительная арифметика — целая наука со своими формулами, аксиомами и теоремами. Целью данной статьи не было давать введение в эту науку. Подходы, приведенные ниже, в первую очередь ориентированы на программистов, решающих инженерные и прикладные задачи, т.е. такие, где диапазоны допустимых значений и необходимые точности вычислений известны и ограничены.
Еще одно ограничение статьи — здесь не приводятся алгоритмы тригонометрических и прочих сложных операций. Сделать их полный обзор в одной статье нереально (и вряд ли необходимо). В статье дан базис, необходимый для понимания таких алгоритмов (и разработки собственных), — правила выполнения базовых операций (сложение/вычитание, умножение, деление) и общая методика вычислений с фиксированной точкой.

Сложение и вычитание

Сложение выполняется просто, если представить в уме, что мы должны сложить две десятичные дроби «в столбик» на листе бумаги. При выполнении данной операции числа записываются в столбик так, чтобы запятые, отделяющие дробную часть, располагались одна под другой. В двоичной арифметике подобная операция называется приведением экспонент.
Если перейти от «бумажки» к математической записи, получится следующее:
Пусть имеется два числа a = n1 * 2^-q1 и b = n2 * 2^-q2.
Тогда:
a + b = n1 * 2^-q1 + n2 * 2^-q2 = (n1 + n2 * 2^{(q1 — q2)})*2^-q1.
Множитель 2^{(q1 — q2)} при втором слагаемом по сути означает арифметический сдвиг для приведения чисел к одной экспоненте.
Стоит отметить, что результат вычисления также может сдвигаться для приведения к требуемому значению экспоненты.
Фрагмент кода на С:

int32_t a = 0x1000L;    // q15: a = 0.125
int32_t b = 0x20000L;   // q20: b = 0.125
int32_t c = 0;          // q25
c = (a << 5) + b;       // q20: (a * 2 ^ (20 - 15) + b); c = 0x40000L (0.25 в q20)
c <<= 5;                // q25: c = 0x800000L (0.25 в q25)

Из примера может быть понятно, что в реальных вычислениях даже в такой простой операции как сложение, есть простор для раздумий. Всегда стоит держать в уме вопросы:

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

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

Умножение

Умножение с фиксированной точкой может выполняться без хитрых выравниваний и приведения к единой экспоненте. Тем не менее умножение — достаточно опасная операция, которая чаще всего в итоге приводит к потере точности и требует особой аккуратности в обращении.
Начнем с математического описания умножения:
Пусть имеется два числа a = n1 * 2^-q1 и b = n2 * 2^-q2.
Тогда:
a * b = n1 * 2^-q1 * n2 * 2^-q2 = n1 * n2 * 2^{-(q2 + q1)}.
Из выражения видно, что экспоненты чисел при умножении складываются: 2^{-(q2 + q1)}. Разрядность данных в этой статье не рассматривается, пока достаточно лишь запомнить, что для безопасного умножения без переполнения и потери точности разрядность результата должна быть не меньше суммарной разрядности сомножителей.
Из-за сложения экспонент результат умножения приходится корректировать для выполнения дальнейших вычислений. При уменьшении экспоненты младшие разряды результата отбрасываются. То есть происходит потеря точности. Можно уменьшить потери точности (и иногда приходится), но способы борьбы с потерями всегда связаны с накладными расходами.
Фрагмент кода на С:

int32_t a = 0x8000L;    // q16: a = 0.5
int32_t b = 0x100000L;  // q21: b = 0.5
int32_t c = 0xC0000L;   // q20: c = 0.75
int64_t d;              // Временная переменная с увеличенным числом разрядов, чтобы хватило на результат.
d = (int64_t)a * (int64_t)b; // q37 = q16 * q21; d = 0x800000000L (0.25 in q37)
d >>= 17;               // q37 / 2 ^ 17 = q20
c += (int32_t)d;        // q20: c = 0x100000 (1 in q20)

Отмечу, что 15 младших разрядов результата умножения были отброшены, чтобы привести число к формату слагаемого. Можно было, конечно, увеличить разрядность переменной c, но, как я уже говорил, на практике диапазоны значений обычно ограничены и младшими разрядами умножения зачастую пренебрегают. Кроме того, не учитывается возможность наличия в исходных сомножителях ненулевых старших разрядов.
Но в данной статье обработка переполнений не рассматривается.

Деление

Начнем с математического выражения для деления:
Пусть имеется два числа a = n1 * 2^-q1 и b = n2 * 2^-q2.
Тогда:
a / b = n1 * 2^-q1 / (n2 * 2^-q2) = n1 / n2 * 2^{-(q1 — q2)}.
Сомножитель 2^{-(q1 — q2)} означает, что при выполнении деления экспонента автоматически уменьшается. Если не принять меры, часть значащих разрядов отбрасывается автоматически.
Способ коррекции очевиден — необходимо заранее увеличить разрядность делителя настолько, чтобы в результате деления получить желаемое количество значащих бит:
a / b = n1 * 2^-q1 * 2^q3 / (n2 * 2^-q2) = n1 / n2 * 2^{-(q1 — q2 + q3)}.
Таким образом, экспонента частного увеличена на q3 разряда.
Фрагмент кода на С:

int32_t a = 0x4000L;    // q15: a = 0.5
int32_t b = 0x80000L;  // q20: b = 0.5
int32_t c = 0;          // q25
int64_t d;              // Временная переменная с увеличенным числом разрядов.
d = (int64_t)a  << 30;  // q45: d = 0x200000000000; (0.5 in q45)
c = (int32_t)(d / (int64_t)b);  // q25: c = 0x2000000; (1 in q25)

Очевидно, что при превышении числом разрядности 32 бита, проблему уже не решить так просто. Тем не менее, для простых инженерных расчетов 32-битных чисел обычно более, чем достаточно.
Есть один простой способ значительно сократить потерю точности при делении — предварительное нормирование делимого. Нормирование — фактически максимальный сдвиг мантиссы влево, при котором не происходит отбрасывания значащих битов. Определить, на сколько можно сдвинуть число, можно путем подсчета ведущих нулей в делимом, для чего существуют специальные алгоритмы (или даже аппаратные инструкции процессора).
После деления частное следует сдвинуть вправо на такое же количество бит для восстановления экспоненты.
Вышеприведенный фрагмент кода при этом может выглядеть таким образом:

int32_t a = 0x4000L;    // q15: a = 0.5
int32_t b = 0x80000L;  // q20: b = 0.5
int32_t c = 0;          // q25
int norm_shift = norm(a); // Вычисление нормирующего сдвига. norm_shift = 16
c = ((a << norm_shift) / b); // q(-5): c = 0x800 (1*2^norm in q(-5))
c <<= (30 - norm);      // q25: c = 0x2000000; (1 in q25)

Как видим, потери точности в данном случае не произошло и без увеличения разрядности делимого.
Впрочем, так происходит не всегда, и если требуется остаться в рамках определенной разрядности (например, 32 бита), приходится реализовывать деление алгоритмически. В обзорной статье вряд ли стоит погружаться в такие дебри — для понимания процесса деления и связанных с ним сложностей достаточно уже приведенного описания.

Нотация, принятая в литературе и различных исходниках

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

QM — где M — число разрядов после запятой. Использована в статье
QN.M — где N — число разрядов до запятой без учета знакового бита, а M — после.

Минус первой нотации очевиден: при работе с переменной приходится обращаться к объявлению переменной (вспоминать ее разрядность) и производить в уме некоторые вычисления, чтобы понять, как привести экспоненту к желаемой. Больше того, если вспомнить округление (int32_t)d в примере с умножением, можно отметить, что при комментариях в данной нотации сложно понять, приведет ли сдвиг или отбрасывание значащих битов к ошибке.
При использовании комментариев во второй нотации можно просто записывать точность вычислений, что исключает необходимость вспоминать, как объявлена переменная.
Поясню примером:

a = 0x1000; // Q15
b = 0x8000; // Q15
int32_t c = a + b; // ??? Без обращения к объявлениям неизвестно, поместится ли результат вычисления в переменную.
a = 0x2000; // Q0.15
b = 0x8000; // Q16.15
c = a + b;  // Q0.15 + Q16.15 = Q16.15: 16 + 15 = 31 бит + 1 знаковый бит

Комментарии во второй нотации очевидно удобнее.
Не буду рассуждать тут о пользе комментариев (далеко не везде они вообще есть), скажу только, что для себя я всегда расписываю при вычислениях типы переменных, чтобы не ошибиться и не запутаться с приведением экспонент.
Если комментарии отсутствуют в принципе, чтение и понимание кода конечно же усложняется, но, запутавшись, всегда можно добавить подобные расшифровки в Q-нотации, чтобы понять, откуда взялся «этот сдвиг влево на 4, а потом сдвиг вправо на 10».
Отмечу, что в недавно выложенных Google исходниках VoIP движка GIPS (webrtc) в комментариях чаще всего просто пишут Q, подразумевая что все биты числа отводятся под дробную часть. Меня лично это весьма путает, т.к. приходится рыться в определениях, чтобы уточнить, как работает код.
Для себя я использую еще одну нотацию, которая отличается от вышеприведенных и близка к нотации MATLAB тулбокса для работы с фиксированной точкой. Она привязывает математику к разрядности переменных и упрощает жизнь, когда надо оценить результат операции (разрядность и экспоненту).
Числа с фиксированной точкой в своих комментариях я отмечаю как QN.M, где N — разрядность числа, M — количество разрядов после запятой.
Поясню, почему я нашел такую схему удобной для себя:

Зная разрядность числа всегда можно предсказать разрядность результата, т.е. выбрать тип переменной, достаточный для его представления.
Мне лично неудобна для чтения запись вида Q(-N).M, которая появляется во второй нотации после выполнения сдвига вправо и нехватке разрядов для дробной части. Например, запись для 16-битного числа, в котором экспонента равна 18 (n*2^-18), у меня выглядит q16.18, а по второй нотации q(-3).18. Запись в первой нотации, как уже было сказано, в любом случае заставляет обращаться к определению для понимания точности вычислений, но в данном случае без определения еще и непонятно: были уже отброшены ведущие значащие биты или нет.
Произведя вычисления по своей нотации, мне проще увидеть, в переменную какой разрядности поместится результат, и как выравнивать экспоненты. Например, q32.15 * q16.4 = q48.19. Сразу видно, что для полного представления результата надо 48-бит. Во второй нотации запись выглядит как q16.15 * q11.4 = q27.19, и приходится подсчитывать, что 27 + 19 = 47 + 1 знаковый от первого сомножителя + 1 знаковый от второго = 48 бит. Мелочь, а приятно. Особенно, когда исходников много.

О плюсах и минусах использования фиксированной точки

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

Необходимость думать.
Предсказуемость результата. При правильном подходе к кодированию результат вычислений будет одинаков на любой платформе (процессор + компилятор) с точностью до разряда. Для данного явления существует специальный термин «битэкзактность» (от англ., bit-exactness). Правильно закодированный алгоритм всегда битэкзактен и, следовательно, может исследоваться на нецелевой платформе. Особенно это полезно, когда отладка на целевой платформе затруднена или невозможна и можно снять только входные данные.
Полный контроль за поведением кода. Фиксированная точка исключает появление «неожиданностей», связанных с особенностями реализации плавающей запятой на используемой платформе.
Автоматическая «фильтрация» пренебрежимо малых значений. В плавающей запятой ошибки вычислений могут накапливаться, в фиксированной точке этого не происходит (за счет отбрасывания малых значений) или процесс накопления ошибок можно контролировать алгоритмически.
Алгоритмически контролируемый диапазон значений переменных. Плавающая запятая дает больше свободы в вычислениях, но результат может выходить за пределы допустимых, что приводит к необходимости его контролировать отдельно. В фиксированной точке эта проблема решается автоматически на этапе разработки и отладки алгоритма.
Переносимость алгоритмов. Данный плюс изрядно коррелирует с первым, но стоит отметить, что целочисленные вычисления гораздо лучше поддержаны множеством не-х86 процессоров, чем вычисления с плавающей запятой. Так что, разработав один раз алгоритм в фиксированной точке, портировать его на различные «слабые» платформы становится гораздо проще. Иногда целочисленные вычисления вообще единственное, что доступно на целевой платформе.
Возможность контролировать сложность вычислений путем понижения точности при разработке алгоритма.
Иногда это интересно.

Минусы:

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

В качестве заключения

Я не стал касаться таких (базовых!) проблем при вычислениях с фиксированной точкой как переполнение результата и дрейф нуля при округлении и методов борьбы с ними. Текст уже получился объемным и, возможно, занудным, чтобы дополнять его подробностями.
В свою очередь я уделил достаточно много внимание общепринятой нотации при записи операций с фиксированной точкой для того, чтобы облегчить чтение специальной литературы и написание собственных исходников (и, возможно, понимание второй части статьи). Да, комментарии с вычислениями в Q-нотации не раз спасали меня от серьезной отладки и разбора исходников.
Если тема будет востребована, я дополню статью следующей частью, в которой опишу вышеуказанные моменты и постараюсь рассказать, как в общем случае можно перевести алгоритм с плавающей точки на фиксированную.
N.B. Математиков просьба не беспокоиться, думаю, вы с освещаемым вопросом знакомы лучше.

Ссылки

UPDATES:

Нотация, которую использую я, как оказалось, пришла мне из MATLAB. Я все не мог вспомнить, откуда я ее выкопал в свое время. Спасибо nerudo, напомнил. Fixed-point toolbox указанного пакета использует для создания объектов как раз пару «число бит» + «число бит на дробную часть», причем с явным указанием: знаковое или беззнаковое число.
Несмотря на то, что примеры изобилуют 8-, 16-, 32- и 64-битными словами, я не старался привязать описание только к х86 и другим General Purpose процессорам. Просто, во-первых так проще приводить примеры на Си, и, во-вторых, я далеко не знаток FPGA, ASIC и пр. (читай VERILOG и т.п.), которые позволяют выбирать произвольную разрядность числа, чтобы внятно о них рассказывать. Возможно, люди, более сведующие в вопросе, смогут дополнить статью своими топиками с примерами.
В описаниях нотаций не сказано о том, что делать с беззнаковыми числами. В общем-то в первоисточниках, в которых я их видел, также ничего не говорится о знаках. В основном подразумевалось, что все числа знаковые. С ходу я не могу сказать, как авторы видели запись знаковых чисел. Отмечу только, что в своей/матлабовской нотации я добавлял букву 'u' после Q, чтобы отметить беззнаковое число.

Tags:

Hubs:

Programming