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

info

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

 

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

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

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

x = m * N e ,

где m и e — целые чис­ла, записан­ные в бинар­ном виде в выделен­ных под них битах. Что­бы избе­жать неод­нознач­ности, счи­тает­ся, что 1 <= |m| < N, то есть чис­ло записа­но в том виде, как если бы оно было с одним зна­ком перед запятой, но запятую злос­тно стер­ли, и чис­ло прев­ратилось в целое.

info

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

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

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

 

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

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

info

Кро­ме оди­нар­ной и двой­ной точ­ности, в новой редак­ции стан­дарта 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, поэто­му в записи ман­тиссы всег­да под­разуме­вает­ся один бит, который не хра­нит­ся в дво­ичной записи. В битах ман­тиссы хра­нит­ся имен­но дроб­ная часть нор­мализо­ван­ного чис­ла в дво­ичной записи.

info

Ман­тисса записы­вает­ся в дво­ичном виде, и отбра­сыва­ется целая часть, заведо­мо рав­ная 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 308 , даже экспо­нен­та оди­нар­ной точ­ности дает воз­можность кодиро­вать чис­ла свы­ше 10 38 . Если срав­нить с «жал­кими» 10 19 , мак­симумом для 64-бит­ных целых чисел, мож­но сде­лать вывод, что мак­сималь­ные и минималь­ные зна­чения вряд ли ког­да‑либо при­дет­ся учи­тывать, хотя и забывать про них не сто­ит.

info

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

Дру­гое дело проб­лема точ­ности. Жал­кие 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 18 .

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

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

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

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

 

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

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

Сло­жение реали­зует­ся срав­нитель­но прос­то:

x = a + b * 10 –18 ,
y = c + d * 10 –18 ,
z = x + y = e + f * 10 –18 ,
a, c, e: int64_t;
b, d ,f: uint64_t;
0 <= b, d, f < 10 18 ,
z = (a + b * 10 –18 ) + (c + d * 10 –18 )
e = a + c + [b * 10 –18 + d * 10 –18 ]
f = {b * 10 –18 + d * 10 –18 } * 10 18

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

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

info

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

Ра­зуме­ется, сто­ит про­верить гра­нич­ные зна­чения при сло­жении a и c. Так­же, исхо­дя из того, что b и d мень­ше 10 18 , мы зна­ем, что (b + d) < 2 * 10 18 , зна­чит, пос­ледним сло­жени­ем добавит­ся мак­симум еди­ница, поэто­му алго­рит­мичес­ки деление мож­но вооб­ще не счи­тать для опти­миза­ции ско­рос­ти выпол­нения опе­раций:

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

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

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

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

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

 

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

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

Ты пом­нишь шко­лу и умно­жение стол­биком? Если нет, самое вре­мя вспом­нить:

a = s a * a 1 - a 2 * 10 -9 ; b = b 1 - b 2 * 10 -9 ;
c = s c * c 1 - c 2 * 10 -9 ; d = d 1 - d 2 * 10 -9 ;
0 <= a 2 , b 2 , c 1,2 , c 1,2 < 10 9 ;
s a,c = sign(a), sign(c)
0 <= a 1 , с 1 < MAX_INT64 / 10 9

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

U = (a 1 , a 2 , b 1 , b 2 ), V = (c 1 , c 2 , d 1 , d 2 ) T , A = V * U,
a 1 *c 1 a 1 *c 1 b 1 *c 1 b 2 *c 1 |
A = | a 1 *c 2 a 1 *c 2 b 1 *c 2 b 2 *c 2 |
a 1 *d 1 a 1 *d 1 b 1 *d 1 b 2 *d 1 |
a 1 *d 2 a 1 *d 2 b 1 *d 2 b 2 *d 2 |

Мат­рица вво­дит­ся не столь­ко для удобс­тва вычис­ления, сколь­ко для удобс­тва про­вер­ки. Ведь A 11 = a 1 * c 1 дол­жно быть стро­го мень­ше MAX_INT64 / 10 18 , а зна­чения диаго­налью ниже: A 12 = a 1 * c 2 и A 21 = a 2 * c 1 дол­жны быть стро­го мень­ше MAX_INT64 / 10 9 . Прос­то потому, что мы будем умно­жать на эти коэф­фици­энты при сло­жении ком­понент:

e = A 11 *10 18 + (A 12 +A 21 )*10 9 + (A 13 +A 22 +A 31 ) + (A 14 +A 23 +A 32 +A 41 ) div 10 18 ,
f = (A 14 +A 23 +A 32 +A 41 ) mod 10 18 + (A 24 +A 33 +A 42 ) + (A 34 +A 43 ) div 10 9

Здесь мы опус­каем сла­гаемое A 44 div 10 18 прос­то потому, что оно рав­но нулю.

Ра­зуме­ется, перед каж­дым сло­жени­ем сто­ит про­верить, не вый­дем ли мы за пре­делы 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 * 10 -18 ) = 10 18 / b,
e = 10 18 div b,
f = 10 18 mod b;

ес­ли b = 0, a = 1, то y = e = 1, f = 0;

ecли b = 0, a > 1, то:
y = 1 / a,
e = 0, f = 10 18 div a.

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

ес­ли a > 1, b != 0 то:
y = 1 / (a + b * 10 -18 ) < 1,

от­сюда e = 0,
f = 10 18 / (a + b * 10 -18 ).

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

k = max(k): 10 k <= a,
u = 10 18 , v = (a * 10 18-k + b div 10 k );
f = (u / v) * 10 18-k ,
for (++k; k <=18; ++k)
{
u = (u % v) * 10;
if (!u) break; // дальше будут нули
f += u / v * 10 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) * 10 18 + 0.5))
{
normalize();
}
void decimal::normalize()
{
uint64_t tail = m_fractional % 10 3 ;
if (tail)
{
if (tail > 10 3 /2)
m_fractional += 10 3 - tail;
else
m_fractional -= tail;
}
}

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

Для float все выг­лядит похожим обра­зом, но пог­решность на порядок выше: 10 18-7 = 10 11. Все целые чис­ла пре­обра­зовы­вают­ся в decimal без проб­лем, прос­то ини­циали­зируя поле m_integral. Пре­обра­зова­ние в обратную сто­рону для целых чисел так­же будет прос­то воз­врат m_integral, мож­но добавить округле­ние m_fractional.

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

return m_integral + m_fractional * 10 -18 ;

От­дель­но сто­ит рас­смот­реть пре­обра­зова­ние в стро­ку и из стро­ки.

Це­лочис­ленная часть, по сути, пре­обра­зует­ся в стро­ку как есть, пос­ле это­го оста­ется толь­ко вста­вить decimal separator и вывес­ти дроб­ную часть как целое, отбро­сив завер­шающие нули.

Так­же мож­но ввес­ти поле «точ­ность» m_precision и записы­вать в стро­ку лишь ука­зан­ное в нем чис­ло десятич­ных зна­ков.

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

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

www

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

 

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

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

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

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

Будь вни­мате­лен к точ­ности, ведь точ­ность навер­няка вни­матель­на к тебе!

Оставить мнение