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

Наглядное объяснение чисел с плавающей запятой

В начале 90-х создание трёхмерного игрового движка означало, что вы заставите машину выполнять почти не свойственные ей задачи. Персональные компьютеры того времени предназначались для запуска текстовых процессоров и электронных таблиц, а не для 3D-вычислений с частотой 70 кадров в секунду. Серьёзным препятствием стало то, что, несмотря на свою мощь, ЦП не имел аппаратного устройства для вычислений с плавающей запятой. У программистов было только АЛУ, перемалывающее целые числа.

При написании книги Game Engine Black Book: Wolfenstein 3D я хотел наглядно показать, насколько велики были проблемы при работе без плавающей запятой. Мои попытки разобраться в числах с плавающей запятой при помощи каноничных статей мозг воспринимал в штыки. Я начал искать другой способ. Что-нибудь, далёкое от и их загадочных экспонент с мантиссами. Может быть, в виде рисунка, потому что их мой мозг воспринимает проще.

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

Как обычно объясняют числа с плавающей запятой

Цитирую Дэвида Голдберта (David Goldbert):

Для многих людей арифметика с плавающей запятой кажется каким-то тайным знанием.

Полностью с ним согласен. Однако важно понимать принципы её работы, чтобы полностью осознать её полезность при программировании 3D-движка. В языке C значения с плавающей запятой — это 32-битные контейнеры, соответствующие стандарту IEEE 754. Они предназначены для хранения и выполнения операций над аппроксимациями вещественных чисел. Пока я видел только такое их объяснение. 32 бита разделены на три части:

Три части числа с плавающей запятой.

Пока всё нормально. Пойдём дальше. Способ интерпретации чисел обычно объясняется с помощью такой формулы:

Именно это объяснение чисел с плавающей запятой все ненавидят.

И здесь я обычно начинаю терять терпение. Возможно, у меня аллергия на математическую нотацию, но когда я это читаю, в моём мозгу ничего не «щёлкает». Такое объяснение похоже на способ рисования совы:

Другой способ объяснения

Хоть это изложение и верно, такой способ объяснения чисел с плавающей запятой обычно не даёт нам никакого понимания. Я виню эту ужасную запись в том, что она разочаровала тысячи программистов, испугала их до такой степени, что они больше никогда не пытались понять, как же на самом деле работают вычисления с плавающей запятой. К счастью, их можно объяснить иначе. Воспринимайте экспоненту как окно (Window) или интервал между двумя соседними целыми степенями двойки. Мантиссу воспринимайте как смещение (Offset) в этом окне.

Три части числа с плавающей запятой.

Окно сообщает нам, между какими двумя последовательными степенями двойки будет число: [0,1], [1,2], [2,4], [4,8] и так далее (вплоть до [,]. Смещение разделяет окно на сегментов. С помощью окна и смещения можно аппроксимировать число. Окно — это отличный механизм защиты от выхода за границы. Достигнув максимума в окне (например, в [2,4]), можно «переплыть» вправо и представить число в пределах следующего окна (например, [4,8]). Ценой этого будет только небольшое снижение точности, потому что окно становится в два раза больше.

Викторина: сколько точности теряется, когда окно закрывает больший интервал? Давайте возьмём пример с окном [0,1], в котором 8388608 смещений накладываются на интервал размером 1, что даёт нам точность . В окне [2048,4096] 8388608 смещений накладываются на интервал , что даёт нам точность .

На рисунке ниже показано, как кодируется число 6,1. Окно должно начинаться с 4 и заканчиваться следующей степенью двойки, т.е. 8. Смещение находится примерно посередине окна.

Значение 6,1 аппроксимированное с помощью числа с плавающей запятой.

Давайте возьмём ещё один пример с подробным вычислением представлением в виде числа с плавающей точкой хорошо известного всем нам значения: 3,14.

Двоичное представление с плавающей точкой числа 3,14.

То есть значение 3,14 аппроксимируется как 3,1400001049041748046875.

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

И, наконец, графическое представление с окном и смещением:

Окно и смещение числа 3,14.

Интересный факт: если модули операций с плавающей запятой были такими медленными, почему в языке C в результате использовали типы float и double? Ведь в машине, на которой изобретался язык (PDP-11), не было модуля операций с плавающей запятой! Дело в том, что производитель (DEC) пообещал Деннису Ритчи и Кену Томпсону, что в следующей модели он будет. Они были любителями астрономии и решили добавить в язык эти два типа.

Интересный факт: те, кому в 1991 году действительно нужен был аппаратный модуль операций с плавающей запятой, могли его купить. Единственными, кому он мог понадобиться в то время, были учёные (по крайней мере, так Intel понимала потребности рынка). На рынке они позиционировались как «математические сопроцессоры». Их производительность была средней, а цена огромной (200 долларов 1993 года — это 350 долларов в 2016 году.). В результате уровень продаж оказался посредственным.

Источник

Разбираемся в числах с плавающей точкой (часть 0)

Здравствуйте, хабровчане. Я давно увлекаюсь темой регистров с плавающей точкой. Меня всегда волновало то, как происходит вывод на экран и т.д. Помню, давным-давно в универе реализовывал свой класс чисел с плавающей точкой, состоящих из 512 бит. Единственное, что я не мог никак реализовать — это вывод на экран.

Как только у меня появилось свободное время, я взялся за старое. Завел себе тетрадку и пошло-поехало. Хотелось додуматься до всего самому, лишь иногда заглядывая в стандарт IEEE 754.
И вот что из всего этого вышло. Интересующихся прошу под кат.

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

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

В сегодняшней статье для примера я буду использовать 32-битные регистры. Числа с двойной точностью (64-битные) работают абсолютно по той же логике.

Сначала поговорим о том, как хранятся числа с плавающей точкой. Старший 31 бит у нас знаковый. Единичка значит, что число отрицательное, а ноль, соответственно наоборот. Далее идут 8 бит экспоненты. Эти 8 битов представляют из себя обычное беззнаковое число. И в самом конце идут 23 бита мантиссы. Для удобства будем обозначать знак как S, экспоненту как E, а мантиссу, как ни странно, M.

Получаем общую формулу

У мантиссы принято считать один неявный единичный бит. То есть мантисса будет представлять из себя 24 бита, но так как старший 23-й бит всегда единица, то его можно не записывать. Таким образом это «ограничение» будет давать нам единственность представления любого числа.

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

Давайте попробуем показать это на наглядном примере:

Представим число 3.625 в двоичном виде. Поначалу разобьем это число на степени двойки.

Степень старшей двойки равна единице. E – 127 = 1. E = 128.

0 1000000 11010000000000000000000

Вот и все наше число.

Попробуем также и в обратную сторону. Пусть у нас есть 32 бита, произвольные 32 бита.

0 10000100 (1)11011100101000000000000

В скобочках указан тот самый неявный старший бит.

Сначала вычислим экспоненту. E = 132. Соответственно степень старшей двойки будет равна 5. Итого имеем следующее число:

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

Для удобной конвертации я накидал небольшую программу на языке C.

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

Можно выразиться иначе. Два соседних числа с плавающей точкой будут отличаться на 2 ^ (E — 127 — 23). То есть на разницу, равную значению младшего бита.

В качестве доказательства можете поменять main в коде и скомпилировать еще раз.

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

Источник

Арифметика с плавающей запятой

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

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

Нужно еще разобраться.

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

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

Что такое мантисса и порядок?

Мантисса – это целое число фиксированной длины, которое представляет старшие разряды действительного числа. Допустим наша мантисса состоит из трех бит (|M|=3). Возьмем, например, число «5», которое в двоичной системе будет равно 101₂. Старший бит соответствует 2 2 =4, средний (который у нас равен нулю) 2 1 =2, а младший 2 0 =1.

Порядок – это степень базы (двойки) старшего разряда. В нашем случае E=2. Такие числа удобно записывать в так называемом стандартном виде, например «1.01e+2». Сразу видно, что мантисса состоит из трех знаков, а порядок равен двум.

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

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

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

Чуть более подробное объяснение:

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

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

Как можно получить NaN? Одним из следующих способов:

Источник

Что нужно знать про арифметику с плавающей запятой

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

Экскурс

Число с плавающей запятой представлено в следующем виде:

$$N = [s] \times m \times 2^,$$

где \(m\) — мантисса ( 23 бита — дробная часть мантиссы и неявно заданный ведущий бит, всегда равный единице, поскольку мантисса хранится в нормализованном виде), \(e\) — смещённая экспонента/порядок ( 8 бит). Один бит отводится под знак ( \(s\), бит равен нулю — число положительное, единице — отрицательное).

Рассмотрим какой-нибудь пример. Число \(3<,>5\) будет представлено в следующем виде:

Бит знака равен нулю, это говорит о том, что число положительное.

Смещённая экспонента равна \(128_<10>\) ( \(10000000_<2>\)), смещена она на \(127\), следовательно, в действительности получаем:

Нормализованная мантисса равна \([1<,>]11000000000000000000000_<2>\), то есть \(1<,>75_<10>\):

$$N = m \times 2^ = 1<,>75 \times 2^ <1>= 3<,>5.$$

Инвертируем бит знака :

Теперь число стало отрицательным: \(-3<,>5\).

Рассмотрим несколько исключений. Как представить число \(0\), если мантисса всегда хранится в нормализованном виде, то есть больше либо равна \(1\)? Решение такое — условимся обозначать \(0\), приравнивая все биты смещённого порядка и мантиссы нулю:

Ноль бывает и отрицательным:

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

В некоторых случаях результат не определён, для этого есть метка \(NaN\), что значит not a number ( не число). Приравняем все биты смещённого порядка и мантиссы единице ( по правде говоря, достаточно того, чтобы хотя бы один бит мантиссы равнялся единице):

Более подробно о типах с плавающей запятой можно прочитать в обучающих статьях. Описанный подход к представлению чисел является компромиссом между точностью, диапазоном значений и скоростью выполнения операций с числами с плавающей запятой. Во-первых, число хранится в двоичной форме, но не всякое десятичное число можно представить в виде непериодической дроби в двоичной системе; это одна из причин потери точности ( \(0<,>1_ <10>= 0<,>000(1100)_<2>\)). Во-вторых, значащих цифр ( цифр мантиссы) не так уж и много: 24 двоичных цифры, или 7,2 десятичных, поэтому придётся много округлять. Но есть и менее очевидные моменты.

Нарушение свойства ассоциативности

Вычислим \(sin(x)\), разложив эту функцию в ряд Тейлора:

Напишем небольшую функцию, которая будет реализовывать эти вычисления. Тип float вместо double выбран для того, чтобы показать, как быстро накапливается погрешность; никаких дополнительных действий не производится, код исключительно демонстрационный. Целью не является показать, что семи членов ряда недостаточно, цель — показать нарушение свойства ассоциативности операций:

хотя свойство коммутативности сохраняется:

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

Посмотрим на это с точки зрения машины:

До этого момента мнение компьютера и человека, казалось бы, сходится. Но результатом в первом случае окажется \(0<,>0000211327\) ( \(2.11327e<->005\)), а во втором случае — \(0<,>0000212193\) ( \(2.12193e<->005\)), совпадают лишь две значащие цифры!

Разгадка проста: у чисел представленного ряда шесть различных ( двоичных) порядков: \(1\), \(2\), \(1\), \(-1\), \(-4\), \(-8\), \(-12\). Когда складываются ( вычитаются) два числа одного порядка или близких порядков, потери точности, как правило, небольшие. Если бы мы складывали огромное число и много маленьких чисел одного порядка, то мы заметили бы, что лучше в плане точности сперва сложить все маленькие числа, а затем уже прибавить большое. Рассмотрим обратный сценарий: сложим большое и первое маленькое число; поскольку порядки значительно различаются, маленькое число будет ( фигурально выражаясь) « раздавлено» большим из-за приведения порядков; получилось новое большое число, не очень точное, но пока ещё достаточно близкое к точному результату; к получившемуся большому числу прибавляем второе маленькое, порядки снова значительно различаются, снова маленькое число оказывается раздавленным, уже две « жертвы». И так далее. Погрешность накопилась достаточно большая.

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

Большие числа

Множество чисел, представимых с помощью float ( и double ), конечно. Из этого, а также из того, что количество разрядов мантиссы весьма ограничено, сразу же следует, что очень большие числа представить не получится; придётся изобретать пути обхода переполнения. Менее очевидным следствием из способа задания числа в виде, описанном выше, является то, что пространство таких чисел не является равномерным. Что это значит?

Возьмём порядок \(e = 0\) и рассмотрим два стоящих подряд числа \(N_<1>\) и \(N_<2>\): у первого мантисса \(m_<1>=[1<,>]00000000000000000000001_<2>\), у второго — \(m_<2>=[1<,>]00000000000000000000010_<2>\).

Конвертер, где можно посмотреть побитовое представление действительного числа

Воспользуемся конвертером: \(N_<1>=1<,>0000001_<10>\), \(N_ <2>= 1<,>0000002_<10>\). Разница между ними равна \(0<,>0000001_<10>\).

Теперь возьмём порядок \(e = 127\) и снова рассмотрим два стоящих подряд числа: у первого мантисса \(m_<3>=[1<,>]00000000000000000000001_<2>\), у второго — \(m_<4>=[1<,>]00000000000000000000010_<2>\).

\(N_<3>=1<,>701412 \times 10^<38>_<10>\), \(N_ <4>= 1<,>7014122 \times 10^<38>_<10>\). Разница между ними равна \(0<,>0000002 \times 10^<38>_<10>\). В рамках этого типа данных нет никакого способа задать некоторое число \(N\), находящееся в интервале между \(N_<3>\) и \(N_<4>\). На секунду, \(0<,>0000002 \times 10^<38>\) — это 20 нониллионов, иначе говоря, двойка и 31 ноль! Маленький шаг для мантиссы и огромный скачок для всего числа.

Каждый из этих диапазонов разбивается на равное количество интервалов. Следовательно, от \(1<,>0\) до \(1<,>9999999\), от \(16<,>0\) до \(31<,>999998\), от \(1<,>7014118 \times 10^<38>\) до \(3<,>4028235 \times 10^<38>\) находится одинаковое количество доступных значений — \((2^ <23>— 2)\) ( поскольку мантисса, за исключением скрытого бита, имеет \(23\) бита; минус сами границы).

Всё сказанное в равной степени касается отрицательных чисел и чисел с отрицательными порядками.

Заключение

Когда в программе используются числа с плавающей запятой, необходимо обращать внимание на операции сложения и вычитания из-за нарушения свойства ассоциативности вследствие ограниченной точности типа float ( и double ). Особенно в том случае, когда порядки чисел сильно различаются. Другой проблемой является большой шаг между ближайшими представимыми числами больши́х порядков. Что и было показано.

Здесь были освещены только эти два аспекта. Также следует помнить о том, что нельзя непосредственно сравнивать два числа с плавающей запятой ( if (x == y) ); что у диапазона есть ограничения; что следует использовать логарифм, когда происходят вычисления с использованием огромных чисел. Ну и совсем не следует забывать о бесконечностях и метке \(NaN\). Из экзотики — флаги компилятора, касающиеся точности промежуточных результатов, из простых вещей ( хотя и не менее коварных, чем что-либо упомянутое выше) — особенности приведения типов. В общем, в выражении « вагон и маленькая тележка» эта заметка — даже не тележка, а маленькая горстка. Маленькая, но достаточно важная.

И последнее. Когда требуется точность вычислений ( например, при финансовых расчётах), следует использовать собственные или сторонние классы чисел с фиксированной запятой. Если речь идёт о языке, поддерживающем такие типы, следует использовать их. Таковым, например, является тип decimal в языке C♯‎.

Источник

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

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

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

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

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

Чтобы избежать неоднозначности, считается, что 1 = 4 503 599 627 370 496 и спокойно вмещает в себя все 32-разрядные целые, давая сбой только на действительно больших 64-разрядных целых (19 десятичных знаков), где погрешность в сотнях единиц уже, как правило, несущественна. Если же нужна большая точность, то мы в данной статье обязательно в этом поможем.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

GITHUB

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

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

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

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

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

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

Источник