Всё, точка, приплыли! Учимся работать с числами с плавающей точкой и разрабатываем альтернативу с фиксированной точностью десятичной дроби

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

Куда уплывает точка

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

Мантисса — это, по сути, число, записанное без точки. Экспонента — это степень, в которую нужно возвести некое число N (как правило, N = 2), чтобы при перемножении на мантиссу получить искомое число (с точностью до разрядности мантиссы). Выглядит это примерно так:

x = m * N ^ e, где m и e — целые числа, записанные в бинарном виде в выделенных под них битах.

Чтобы избежать неоднозначности, считается, что 1

Мантисса — это, по сути, число, записанное без точки. Экспонента — это степень, в которую нужно возвести некое число N (как правило, N = 2), чтобы при перемножении на мантиссу получить искомое число

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

Современные вещественные числа, поддержанные аппаратно на уровне процессора, также разбиты на мантиссу и экспоненту. Разумеется, все операции, привычные для арифметики целых чисел, также поддержаны командами процессора для вещественных чисел и выполняются максимально быстро.

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

Стандарт точечного плаванья

Вещественные числа с плавающей точкой, поддержанные на уровне процессора, описаны специальным международным стандартом IEEE 754. Основными двумя типами для любых вычислений являются single-precision (одинарной точности) и double-precision (двойной точности) floating-point (числа с плавающей точкой). Названия эти напрямую отражают разрядность бинарного представления чисел одинарной и двойной точности: под представление с одинарной точностью выделено 32 бита, а под двойную, как ни странно, 64 бита — ровно вдвое больше.

Кроме одинарной и двойной точности, в новой редакции стандарта IEEE 754—2008 предусмотрены также типы расширенной точности, четверной и даже половинной точности. Однако в C/C++, кроме float и double, есть разве что еще тип long double, упорно не поддерживаемый компанией Microsoft, которая в Visual C++ подставляет вместо него обычный double. Ядром процессора в настоящий момент также, как правило, пока не поддерживаются типы половинной и четверной точности. Поэтому, выбирая представления с плавающей точкой, приходится выбирать лишь из float и double.

В качестве основания для экспоненты числа по стандарту берется 2, соответственно, приведенная выше формула сводится к следующей:

x = m * (2 ^ e), где 1 <= |m| < 2; m, e — целые

Расклад в битах в числах одинарной точности выглядит так:

| 1 бит под знак | 8 бит экспоненты | 23 бита мантиссы |

Для двойной точности мы можем использовать больше битов:

| 1 бит под знак | 11 бит экспоненты | 52 бита мантиссы |

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

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

Мантисса записывается в двоичном виде, и отбрасывается целая часть, заведомо равная 1, поэтому никогда не забываем, что мантисса на один бит длиннее, чем в она хранится в двоичном виде

Не нужно иметь докторскую степень, чтобы вычислить точность в десятичных знаках чисел, которые можно представить этим стандартом: 2^(23 + 1) = 16 777 216; это явно указывает нам на тот факт, что точность представления вещественных чисел с одинарной точностью достигает чуть более семи десятичных знаков. Это значит, что мы не сможем сохранить в данном формате, например, число 123 456,78 — небольшое, в общем-то, число, но уже начиная с сотой доли мы получим не то число, что хотели. Ситуация усложняется тем, что для больших чисел вида 1 234 567 890, которое прекрасно помещается даже в 32-разрядное целое, мы получим погрешность уже в сотнях единиц! Поэтому, например, в C++ для вещественных чисел по умолчанию используется тип double. Мантисса числа с двойной точностью уже превышает 15 знаков: 2^52> = 4 503 599 627 370 496 и спокойно вмещает в себя все 32-разрядные целые, давая сбой только на действительно больших 64-разрядных целых (19 десятичных знаков), где погрешность в сотнях единиц уже, как правило, несущественна. Если же нужна большая точность, то мы в данной статье обязательно в этом поможем.

Теперь что касается экспоненты. Это обычное бинарное представление целого числа, в которое нужно возвести 10, чтобы при перемножении на мантиссу в нормализованном виде получить исходное число. Вот только в стандарте вдобавок ввели смещение, которое нужно вычитать из бинарного представления, чтобы получить искомую степень десятки (так называемая biased exponent — смещенная экспонента). Экспонента смещается для упрощения операции сравнения, то есть для одинарной точности берется значение 127, а для двойной 1023. Все это звучит крайне сложно, поэтому многие пропускают главу о типе с плавающей точкой. А зря!

Примерное плаванье

Чтобы стало чуточку понятнее, рассмотрим пример. Закодируем число 640 (= 512 + 128) в бинарном виде как вещественное число одинарной точности:

• число положительное — бит знака будет равен 0;
• чтобы получить нормализованную мантиссу, нам нужно поделить число на 512 — максимальную степень двойки, меньшую числа, получим 640 / 512 = 512 / 512 + 128 / 512 или 1 + 1/4, что дает в двоичной записи 1,01, соответственно, в битах мантисы будет 0100000 00000000 00000000;
• чтобы получить из 1 + 1/4 снова 640, нам нужно указать экспоненту, равную 9, как раз 2^9 = 512, число, на которое мы поделили число при нормализации мантиссы, но в бинарном виде должно быть представление в смещенном виде, и для вещественных чисел одинарной точности нужно прибавить 127, получим 127 + 9 = 128 + 8, что в бинарном виде будет записано так: 10001000.

Для двойной точности будет почти все то же самое, но мантисса будет содержать еще больше нулей справа в дробной части, а экспонента будет 1023 + 9 = 1024 + 8, то есть чуть больше нулей между старшим битом и числом экспоненты: 100 00001000. В общем, все не так страшно, если аккуратно разобраться.

Задание на дом: разобраться в двоичной записи следующих констант: плюс и минус бесконечность (INF — бесконечность), ноль, минус ноль и число-не-число (NaN — not-a-number).

За буйки не заплывай!

Есть одно важное правило: у каждого формата представления числа есть свои пределы, за которые заплывать нельзя. Причем обеспечивать невыход за эти пределы приходится самому программисту, ведь поведение программы на С/С++ — это сделать невозмутимое лицо при выдаче в качестве сложения двух больших положительных целых чисел маленькое отрицательное. Но если для целых чисел нужно учитывать только максимальное и минимальное значение, то для вещественных чисел в представлении с плавающей точкой следует больше внимания обращать не столько на максимальные значения, сколько на разрядность числа. Благодаря экспоненте максимальное число для представления с плавающей точкой при двойной точности превышает 10³⁰⁸, даже экспонента одинарной точности дает возможность кодировать числа свыше 10³⁸. Если сравнить с «жалкими» 10¹⁹, максимумом для 64-битных целых чисел, можно сделать вывод, что максимальные и минимальные значения вряд ли когда-либо придется учитывать, хотя и забывать про них не стоит.

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

Другое дело проблема точности. Жалкие 23 бита под мантиссу дают погрешность уже на 8-м знаке после запятой. Для чисел с двойной точностью ситуация не столь плачевная, но и 15 десятичных знаков очень быстро превращаются в проблему, если, например, при обработке данных требуется 6 фиксированных знаков после точки, а числа до точки достаточно большие, под них остается всего лишь 9 знаков. Соответственно, любые многомиллиардные суммы будут давать значительную погрешность в дробной части. При большой интенсивности обработки таких чисел могут пропадать миллиарды евро, просто потому, что они «не поместились», а погрешность дробной части суммировалась и накопила огромный остаток неучтенных данных.

Если бы это была только теория! На практике не должно пропадать даже тысячной доли цента, погрешность всех операций должна быть строго равна нулю. Поэтому для бизнес-логики, как правило, не используют C/C++, а берут C# или Python, где в стандартной библиотеке уже встроен тип Decimal, обрабатывающий десятичные дроби с нулевой погрешностью при указанной точности в десятичных знаках после запятой. Что же делать нам, программистам на C++, если перед нами стоит задача обработать числа очень большой разрядности, при этом не используя высокоуровневые языки программирования? Да то же, что и обычно: заполнить пробел, создав один небольшой тип данных для работы с десятичными дробями высокой точности, аналогичный типам Decimal высокоуровневых библиотек.

Добавим плавающей точке цемента

Пора зафиксировать плавающую точку. Поскольку мы решили избавиться от типа с плавающей точкой из-за проблем с точностью вычислений, нам остаются целочисленные типы, а поскольку нам нужна максимальная разрядность, то и целые нам нужны максимальной разрядности в 64 бита.

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

Отсыпь и мне децимала!

Сначала немного теории. Обозначим наше две компоненты, целую и дробную часть числа, как n и f, а само число будет представимо в виде

x = n + f * 10^-18, где n, f — целые, 0 <= f < 10¹⁸.

Для целой части лучше всего подойдет знаковый тип 64-битного целого, а для дробной — беззнаковый, это упростит многие операции в дальнейшем.

class decimal
{
        …
private:
        int64_t  m_integral;
        uint64_t m_fractional;
};

Целая часть в данном случае — максимальное целое, меньшее представляемого числа, дробная часть — результат вычитания из этого числа его целой части, помноженной на 10¹⁸, и приведенное к целому: f = (x – n) * 10¹⁸.

Целая часть для отрицательных чисел получится большей по модулю самого числа, а дробная часть будет не совпадать с десятичной записью самого числа, например для числа –1,67 компонентами будут: n = –2 и f = 0,33 * 10¹⁸. Зато такая запись позволяет упростить и ускорить алгоритмы сложения и умножения, поскольку не нужно ветвления для отрицательных чисел.

Операции с типом десятичной дроби

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

x = a + b * 1e-18,
y = c + d * 1e-18,
z = x + y = e + f * 1e-18,

a, c, e: int64_t;
b, d ,f: uint64_t;
0 <= b, d, f < 1e+18,

z = (a + b * 1e-18) + (c + d * 1e-18)
e = a + c + [b * 1e-18 + d * 1e-18]
f = {b * 1e-18 + d * 1e-18} * 1e+18

NB: здесь и далее все записи в форме 1e — целые числа.

Здесь [n] — это получение целой части числа, а {n} — получение дробной части. Все бы хорошо, но вспоминаем про ограничение целых чисел. Значение 1e+18 уже близко к грани значений беззнакового 64-битового целого типа uint64_t (потому мы его и выбрали), но нам никто не мешает чуточку упростить выражение, чтобы гарантированно оставаться в границах типа, исходя из начальных условий:

e =  a + c + (b + d) div 1e+18,
f = (b + d) mod 1e+18.

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

Разумеется, стоит проверить граничные значения при сложении a и c. Также, исходя из того, что b и d меньше 1e+18, мы знаем, что (b + d)

e = a + c;
f = b + d;
if (f >= 1e+18) f -= 1e+18, ++e;

Здесь опущены проверки на максимальное целое для значения e для простоты изложения.

Для вычитания все чуть-чуть сложнее, здесь уже нужно учитывать переход ниже нуля для беззнакового целого. То есть нужно сравнить два числа до вычитания.

e = a - c;
if (b >= d) f = b - d;
else f = (1e+18 - d) + b, --e;

В целом все пока несложно. До умножения и деления все всегда несложно.

Умножение чисел с фиксированной точностью

Рассмотрим сперва умножение. Поскольку числа в дробной части довольно большие и, как правило, находятся в ближайшей окрестности 1e+18, нам придется либо использовать язык ассемблера и операции с Q-регистрами, либо пока обойтись целыми числами, побив каждую компоненту на две части по 1e+9. В этом случае умножение будет давать не больше 1e+18 и ассемблер нам пока не понадобится, будет не так шустро, но нам хватит 64-разрядного целого, и мы останемся внутри C++. Ты помнишь школу и умножение столбиком? Если нет, самое время вспомнить:

a = s_a * a₁ — a₂ * 10^-9; b = b₁ — b₂ * 10^-9;
c = s_c * c₁ — c₂ * 10^-9; d = d₁ — d₂ * 10^-9;
0 <= a₂, b₂, c_1,2, c_1,2 < 10⁹;
s_a,c = sign( a ), sign( c )
0 <= a₁, с₁ < MAX_INT64 / 10⁹

Введем матрицу для упрощения вычисления умножения:

U = (a₁, a₂, b₁, b₂),
V = (c₁, c₂, d₁, d₂)^T,
A = V * U,
A=
| a₁*c₁ a₁*c₁ b₁*c₁ b₂*c₁ |
| a₁*c₂ a₁*c₂ b₁*c₂ b₂*c₂ |
| a₁*d₁ a₁*d₁ b₁*d₁ b₂*d₁ |
| a₁*d₂ a₁*d₂ b₁*d₂ b₂*d₂ |

Матрица вводится не столько для удобства вычисления, сколько для удобства проверки. Ведь A₁₁ = a₁ * c₁ должно быть строго меньше MAX_INT64 / 10¹⁸, а значения диагональю ниже: A₁₂ = a₁ * c₂ и A₂₁ = a₂ * c₁ должны быть строго меньше MAX_INT64 / 109. Просто потому, что мы будем умножать на эти коэффициэнты при сложении компонент:

e = A₁₁*10¹⁸ + (A₁₂+A₂₁)*10⁹ + (A₁₃+A₂₂+A₃₁) + (A₁₄+A₂₃+A₃₂+A₄₁) div 10¹⁸,
f = (A₁₄+A₂₃+A₃₂+A₄₁) mod 10¹⁸ + (A₂₄+A₃₃+A₄₂) + (A₃₄+A₄₃) div 10⁹

Здесь мы опускаем слагаемое A₄₄ div 10¹⁸ просто потому, что оно равно нулю. Разумеется, перед каждым сложением стоит проверить, не выйдем ли мы за пределы MAX_INT64. К счастью, мы можем оперировать беззнаковым типом uint64_t для всех компонент матрицы и для промежуточного результата. Все, что нужно будет сделать в конце, — это определить знак результата s_e = s_a xor s_c и для отрицательного числа поправить целую и дробную часть: целую уменьшить на единицу, дробную вычесть из единицы. Вот, в общем, и все умножение, главное — быть очень аккуратным. С ассемблером все на порядок проще, но этот материал выходит за рамки Академии C++.

Алгоритм деления без регистрации и СМС

Если ты помнишь алгоритм деления столбиком — молодец, но здесь он будет не нужен. Благодаря математике и небольшому колдовству с неравенствами нам будет проще посчитать обратное число x^–1 и выполнить умножение на x^–1. Итак, решаем задачу

y = x^-1 = 1 / (a + b * 10^-18) = c + d * 10^-18.

Для упрощения рассмотрим нахождение обратного числа для положительного x. Если хотя бы одна из компонент x равна нулю (но не обе сразу), вычисления сильно упрощаются. Если a = 0, то:

y = 1 / (b * 1e-18) = 1e+18 / b,
e = 1e+18 div b,
f = 1e+18 mod b;
если b = 0, a = 1, то y = e = 1, f = 0;
ecли b = 0, a > 1, то:
y = 1 / a,
e = 0, f = 1e+18 div a.

Для более общего случая, когда x содержит ненулевые дробную и целую части, в этом случае уравнение сводится к следующему:

если a > 1, b != 0 то:
y = 1 / (a + b * 1e-18) < 1, 
отсюда e = 0,
f = 1e+18 / (a + b * 1e-18).

Теперь нужно найти максимальную степень 10, которая будет не больше a, и итерационно выполнять следующее действие:

k = max(k): 1e+k <= a,
u = 1e+18, v = (a * 1e+18-k + b div 1e+k);
f = (u / v) * 1e+18-k,
for (++k; k <=18; ++k)
{
        u = (u % v) * 10;
        if (!u) break; // дальше будут нули
        f += u / v * 1e+(18-k);
}

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

Очень полезно будет завести массив степеней десяток от 0 до 18 включительно, поскольку вычислять их совершенно излишне, мы их знаем заранее и требоваться они нам будут часто.

Преобразования типов

Мы знаем и умеем достаточно, чтобы теперь превратить расплывчатые float и double в наш новенький decimal.

decimal::decimal(double value)
   :    m_integral(static_cast<int64_t>(std::floor(value))
        m_fractional(static_cast<int64_t>(std::floor(
            (value - m_integral) * 1e+18 + 0.5))
{
    normalize();
}

void decimal::normalize()
{
    uint64_t tail = m_fractional % 1e+3;
    if (tail)
    {
        if (tail > 103/2)
            m_fractional += 1e+3 - tail;
        else
            m_fractional -= tail;
    }
}

Здесь 103 является, по сути, той погрешностью, за которой double перестает быть точным. При желании погрешность можно еще уменьшить, здесь 10^18-15 нужно для наглядности изложения. Нормализация после преобразования нужна будет все равно, поскольку точно double заведомо ниже даже дробной части decimal. Кроме того, нужно учитывать случай, когда double выходит за пределы int64_t, при таких условиях наш decimal не сможет правильно преобразовать целую часть числа.

Для float все выглядит похожим образом, но погрешность на порядок выше: 10^18-7 = 10¹¹.

Все целые числа преобразовываются в decimal без проблем, просто инициализируя поле m_integral. Преобразование в обратную сторону для целых чисел также будет просто возврат m_integral, можно добавить округление m_fractional.

Преобразование из decimal в double и float сводится к вышеуказанной формуле:

return m_integral + m_fractional * 1e-18;

Отдельно стоит рассмотреть преобразование в строку и из строки. Целочисленная часть, по сути, преобразуется в строку как есть, после этого остается только вставить decimal separator и вывести дробную часть как целое, отбросив завершающие нули. Также можно ввести поле «точность» m_precision и записывать в строку лишь указанное в нем число десятичных знаков.

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

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

GITHUB

Со статьей идет несколько файлов с исходниками одной из возможных реализаций decimal, а также с небольшим тестом вещественных чисел для лучшего усвоения материала.

Не уплывай, и точка!

В заключение скажу лишь то, что подобный тип в C/C++ может появиться в весьма специфической задаче. Как правило, проблемы чисел с большой точностью решаются языками типа Python или C#, но если уж понадобилось по 15–18 знаков до запятой и после, то смело используй данный тип.

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

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

Впервые опубликовано в журнале Хакер #192.
Автор: Владимир Qualab Керимов, ведущий С++ разработчик компании Parallels

Подпишись на «Хакер»

This entry passed through the Full-Text RSS service - if this is your content and you're reading it on someone else's site, please read the FAQ at http://ift.tt/jcXqJW.