Плаваюча (рухома) крапка, частина 3: Сучасний стандарт і все навколо — основи
Читайте також:
— Плаваюча (рухома) крапка. Частини
— Плаваюча (рухома) крапка, частина 4: Сучасний стандарт — тонкі деталі і проблеми
Тут і далі, знак «↑» позначує піднесення у степінь; так, 2↑5 = 32, e↑(i×π) = —1.
Сучасний стандарт і все навколо
Мабуть, кожен, хто використовував рухому крапку не тільки для цілих чисел, десь чув, що всі сучасні реалізації рухомої крапки відповідають IEEE754. Більш досвідчені можуть ще згадати версії стандарту, поступове введення окремих рис, консенсус між реалізаціями, як він досягався... Складно зараз найти реалізацію архітектури компʼютера, що б не відповідала, згідно твердженням її розробників чи власників, IEEE754. Нових таких не виникає десь з 1990 року, чи ще раніше. (Навіть суперкомпактні формати для даних з ML або нейромереж задовільняють тим же принципам.) Що тут ще казати?
Два моменти чомусь дуже часто втрачаються з уваги, а вони тут принципові.
Перший: IEEE754 — це стандарт IEEE. Не ISO, не ANSI, не IETF, а IEEE — орґанізації, яка має головним націленням практики електроніки, тобто апаратного забезпечення, «хардвера». IEEE має авторитетні результати і в «софтвері», але коріння дають конкретний підхід до всього.
(Зовсім «юристи» можуть згадати, що він же прийнятий як ISO/IEC 60559. Але це тільки «імпорт» готового. Ні ISO ані IEC серйозно руку до нього не прикладали.)
IEEE754 розроблений під те, щоб надати універсальний спосіб зробити оптимальну апаратну реалізацію — в тих умовах, в яких взагалі апаратна реалізація має сенс, і до тих меж, що вона здатна забезпечити.
(Я б сказав «цілком і повністю» розроблений... але це вже не так. Оптимізація апаратної реалізації розрахована на найбільш прості, базові операції — тобто чотири арифметичних, і допоміжні операції. Це вже не розповсюджується на конверсії між внутрішнім двійковим і зовнішнім десятковим представленням, на трансцедентні функції... Але все одно вони базуються на тому ж підґрунті.)
Другий момент: відповідність може бути різна. Фактично, можна сказати, що IEEE754 це не фіксований стан, а процес. І повної відповідности, на диво, нема ні у кого. Або я таких не знаю, біля мене взагалі не пробігали. Але повної відповідности і не треба, якщо підходити з практичних потреб.
Ніхто не зобовʼязаний завжди використовувати тільки те, що є за цим стандартом. Кожен є вільним відмовитись від нього і використовувати щось інше. Наприклад, у embedded типові ситуації, що такий sin(), який в стандарті, просто не потрібний, замість цього треба 255*sin(x*pi/11520) у 200 значень з дозволеною похибкою до 6.5% у табличному вигляді (бо нема часу рахувати у рантаймі), яка дозволяє стиснути представлення в памʼяті. Або незалежна реалізація реґулярно використовується, хоча б для контролю штатних реалізацій. Якщо ви не впевнені в якости того ж sin() чи навіть множення у вашому рантаймі, беріть, наприклад, MPFR, задавайте потрібну точність і звіряйте результати з тим, що дає бібліотека.
Але... саморобні, суто проґрамні аналоги повної функціональности і точности (базових операцій; це вже не завжди для складних функцій) будуть значно повільнішими. Разів в
Ефективна апаратна реалізація в наш час це лише витрата не надмірної кількости грошей. Але критично те, що стандарт дає засоби отримувати однакові результати на різних мовах, компіляторах, процесорах (одної архітектури чи різних). В умовах, коли не можеш зрозуміти, хто винен, що у тебе виник інший розвʼязок рівняння чи обчислення перестали сходитись до якогось конкретного результату, можливість відтворення стає критичною. Зрозуміло, що цими засобами ще треба скористатись — і правильно скористатись.
Після такого вступу можна вже перейти до обговорення деталей стандарту.
В давні часи в староглиняному царстві...
Що було до прийняття IEEE754? Сказати, що був бардак — нічого не сказати. Ось ця сторінка показує тільки формати представлення чисел, яких більше 50! Але вже тут багато цікавого, як можна було майже кожний аспект зробити дещо інакше.
А ще є локальні правила виконання операцій, які теж були свої у кожної моделі компʼютера, і добре якщо не у кожної субмодифікації... Далі буде дещо більше прикладів, але некоректне (або повністю відсутнє) округління, немонотонні результати, нестійка похибка, яка приводила до нестійкости обчислень — все це було проблемами, через які кожна міґрація між моделями компʼютерів була маршем через нерозмічене мінне поле з непередбачуваним часом на стабілізацію після переходу.
Розуміння, що «так жити далі не можна», прийшло фінально в середині
But the cumulative effect of those idiosyncrasies and the programming contortions they induce imposes a numbing intellectual burden upon the software industry.
<...>
These developments could stimulate various imaginings: that computer arithmetic had been in a state of anarchy; that the production and distribution of portable numerical software had been paralyzed that numerical analysts had been waiting for a light to guide them out of chaos.
Почитайте решту, кому цікаво — вся стаття характерна як «дух часу» і дає детальні розповіді деяких дуже специфічних аспектів (котрі я не описуватиму тут через складність).
Драфт стандарту IEEE вже існував, але до прийняття навіть такої скороченої версії, як 1985, було ще 4 роки. Що таке виробити стандарт там, де найменша зміна це витрати на міліони у кожного вендора, може розповісти той, хто читає потік робочих документів якогось комітету з відкритих, як WG21 (C++).
Прийняття IEEE754 було кроком відчаю в цьому всьому. І добре, що змогли прийняти те, що було розумно в основі, зуміло вижити і розвиватись.
Ми ще згадуватимо ці деталі.
Вимоги, версії, варіації стандарту
Є версії 1985, 2008 і 2019 року. 2019 це невеличкі правки до 2008, їх далі розглядатимемо разом, крім випадків, коли ці правки будуть важливі. А ось між 1985 і 2008 — прірва. Версія 1985 року була дуже мінімальною: тільки зовнішнє представлення, найпростіша арифметика, округлення, контекст і помилки; усталена практика вже була значно ширша, ніж стандарт. Навпаки, версія 2008 року (і її розвиток 2019 року) визначає майже все, навіть те, чого не роблять. Орієнтуємось на версію 2019 року, але маємо на увазі все інше.
Крім того, є побічний стандарт IEEE854-1987, який задає деякі загальні характеристики одночасно для основ 2 і 10, але без фіксованих значень. Він цікавий тим, що задає деякі обмеження на чисельні значення, які показують і походження чисельних значень у головному стандарті... але без розʼяснень. Це спільна проблема обох стандартів — відсутність пояснень вибраних рішень: їх треба зʼясовувати з супутних документів і навіть мемуарів розробників. IEEE854 формально скасований, бо увійшов в IEEE754-2008... але це не зовсім правда, якщо прийняти до уваги ті норми IEEE854, які загальніші, ніж треба для 754; він продовжує мати суттєве історичне значення і ми його згадуватимо.
Окремо є набір стандартів language-independent arithmetic (LIA), ISO/IEC 10967. Вони задають рекомендації щодо як арифметики з плаваючою крапкою, так і з цілими числами. Подробний розбір цих стандартів за межами даного циклу, обмежимось окремими згадками. І ще є General Decimal Arithmetic Specification від IBM, про неї буде в частині про десяткову арифметику.
Що визначає стандарт (тут і далі, IEEE754-2019, якщо не уточнено):
1. Зовнішні представлення чисел (interchange formats). Вони не обовʼязкові для внутрішнього використання, але на практиці майже всі використовують саме їх. Є версії для основи 2 і для основи 10; останні зʼявились тільки з версії 2008, мають, на жаль, дуже мале поширення.
2. Розширені представлення/формати (extended і extendable). Не наше питання зараз, але цікавим є співвідношення до відповідних базових форматів.
3. Головні принципи виконання операцій: контекст (режими, винятки), округлення, інші деталі.
4. Обовʼязковий набір операцій. Це по більшости арифметика, округлення, конверсії внутрішніх типів (між плавучими форматами, а також в ціле і з цілого), і конверсії між внутрішнім двійковим і зовнішнім десятковим представленнями.
5. Рекомендований набір операцій (на додачу до обовʼязкового). Це по більшости степені-лоґарифми і триґонометрія (хтось з тих, хто це читає, використовував гіперболічний косінус?)
6. Обробка помилок.
7. Також різне допоміжне (згадаємо, якщо буде сенс).
Фактично, це той же список, що був у першій частині, але переґрупований і з іншими деталями.
Що незручно для аналіза, це те, що в стандарті майже нема обґрунтувань вибору рішення, константи якогось обмеження, інших суттєвих рішень. Інтернет заповнений питаннями «чому саме так?» На щастя, на більшість (не всі) є відповіді, але їх треба розуміти в контексті конкретного часу.
Пройдемось поступово по перелічених аспектах, згадуючи не тільки стандартні вимоги і правила, а ще й якесь обґрунтування (ну... у межах, щоб не переповнити). Але перед цим невеличкий відступ.
Куди IEEE754 не дотягнувся?
Є величезна і вкрай впливова область, де цей стандарт не прийняли і, схоже, ще довго не прийматимуть: SQL і взагалі реляційні бази даних. Там все погано, але про це окрема розмова.
Фінансова арифметика виходить з принципово іншого підходу (в частині 1 вже було: похибка заборонена крім спеціально відведених і позначених місць), і не приймає користь від основи 2. Стандартні норми десяткової арифметики можуть бути застосовані, якщо є. Про неї також буде окрема розповідь.
Представлення чисел
Зараз обмежуємось форматами з основою 2.
Всі транспортні (interchange) формати представлення за IEEE754 мають спільний принцип, відрізняючись тільки константами розмірів: 1 біт на знак, декілька на зсунутий (зміщений) порядок (biased exponent), і решта на significand (еквівалент мантиси).
Базові стандартизовані розміри:
|Формат|Be|Bsg| Мін.дод. |Мін.дод.норм.| Макс.дод. |min10d|max10d| |------|--|---|-------------|-------------|-------------|------|------| |half | 5| 10|5.9605e-08 |6.1035e-05 |65504 | 3 | 5 | |single| 8| 23|1.4013e-45 |1.17549e-38 |3.40282e+38 | 6 | 9 | |double|11| 52|4.94066e-324 |2.22507e-308 |1.79769e+308 | 15 | 17 | |quad |15|112|6.47518e-4966|3.3621e-4932 |1.18973e+4932| 33 | 36 |
(точність значень урізана, інакше од цифр миготить в очах)
Be — кількість бітів порядка; Bsg — significandʼа (мантиса нормалізованого числа має на 1 біт більше, відповідно, 11, 24, 53, 113); min10d — кількість десяткових цифр, ґарантовано збережених у двійковому форматі; max10d — кількість десяткових цифр, потрібних для ґарантованого збереження значення у двійковому форматі. Розшифровка скорочень: мінімальне додатнє (субнормальне); мінімальне додатнє нормалізовано; максимальне додатнє (не нескінченність).
У зсунутого порядку два спеціальні значення: 1) всі біти 0 — на субнормальні числа і точний нуль, і 2) всі біти 1 — на нескінченність (infinity) і не-число (not a number, NaN). Всі інші значення (тобто не крайні) вказують «звичайне», нормалізоване число. Наприклад, для одинарної точности (32 біти повного представлення) відведено 8 біт на поле зсунутого порядку (від 0 до 255), і тоді значення 1...254 вказує на нормалізоване число, 0 — на нуль або субнормальне, 255 — на нескінченність або не-число. Що за «не-число» і «нескінченність», звідки, що і чому з ними робити — буде далі в цій частині і наступних.
Що лежить у significand? Тут правило таке:
1) Якщо порядок показує нормалізоване число, то мантиса нормалізована (1 ⩽ M < 2), старший біт не зберігається (завжди рівний 1), в significand лежить дробова частина мантиси в основі 2.
Наприклад, significand, рівний 00...0 (всі біти 0), показує, таким чином, число, що є степенем двох. Якщо він дорівнює 10...00 — то мантиса дорівнює 1.5 (бо 1.10...00₂ == 1.1₂ == 1.5₁₀)... В загальному випадку, поставте «1.» перед ним і переведіть з двійкової системи.
Його ж можна переінтерпретувати (подумки змістивши порядок) як дріб (звичайний, двох цілих чисел). Наприклад, 0.1 в подвійній точності це 7 205 759 403 792 794 / 72 057 594 037 927 936. Знаменник тут дорівнює 2↑56, а числівник це всі біти significand з доданою 1 спереди. Якщо експортувати в ASN.1 або CBOR (навмисно в bigfloat), це варіант того, що ви можете отримати в двійковому вигляді; там представлення двійкового плаваючого числа це мантиса як ціле число, і до неї порядок; для 0.1 взятого з IEEE double буде мантиса рівна числівнику цього дробу, і порядок —56.
2) Якщо зсунутий порядок дорівнює 0, то significand сприймається як дробова частина при цілій частині мантиси 0, а зсунутий порядок як був би рівний 1. Якщо significand не 0, то це «субнормальне» число. Якщо significand дорівнює 0, то число — просто нуль. В стандарті зараз це зветься subnormal number, раніше був також термін denormalized number, але він комусь здався менш прийнятним. Субнормальні мають скорочену точність, в залежности від їх значення (де найстарша одиниця).
Про субнормальні числа буде далі. Взагалі про них багато міфів. Поки що просто знаємо, що це варіант норми:)
3) Якщо зсунутий порядок — всі біти 1, то це або нескінченність або не-число.
Що значить порядок зсунутий? Зміщення треба відняти зі зсунутого порядку, щоб отримати незсунутий, який показує реальне значення в базовій формулі. Наприклад, зміщення 127 означає, що якщо в полі порядку записано 131, то незсунутий відповідно дорівнює 4, і мантису, отриману з significandʼу (1 ⩽ M < 2), треба домножити на 2↑4. Чому порядок зсунутий? Ця манера виникла ще у перші роки розробки такого апаратного забезпечення. Простіше працювати з числами без знаку, ніж зі знаком — коли вони не виходять за межі.
Але ось чому дорівнює зсув, це значно цікавіше і демонструє деякі глибинні деталі і задуми стандартизаторів... Якщо подивитись в приведений список форматів, типовий зсув (bias в стандарті, excess в списку з quadibloc), яке додано до реального порядку при збереженні у памʼяті, у більшости дорівнює рівно половині діапазона значень (при 8 бітах на зсунутий порядок це 128). Так робити дешево, не було причини робити інакше, і так робили майже всі. Але для IEEE754 приведено 126. Що не так? Це цікаве питання — і воно показує лоґіку стандартизаторів і якість адаптації до «усереднених» потреб.
В першій частині вже вводили деякі потрібні тут терміни. Нехай наше число 12340. Якщо його записати як 0.1234×10↑5 (1/b ⩽ M < 1, де b — основа), це fraction view. 1.234×10↑4 (1 ⩽ M < b) — left units view. Нарешті, 1234×10↑1 чи 12340×10↑0... (M — ціле число; вище ще був приклад для 0.1) — right units view. Всі три форми мають своє призначення. Зараз ми говоримо про left units view і відповідну інтерпретацію незсунутого порядку. У тому, що описується в цьому розділі, це left units view, в ньому простіше за все поясняти.
Зсув має бути максимально близьким до половини діапазону порядків, щоб для якомога більшої кількости чисел одночасно x і 1/x були представлені. Але при розробці IEEE754 додали вимогу: зворотнє до мінімального (по модулю) нормального числа має бути представленим як скінченне число. При цьому дозволено, щоб зворотнє до максимального нормального було вже субнормальним, але таким, що має представлення. Легка втрата точности тут менше зло, ніж вихід у нескінченність, після якого повернутись назад до скінченних значень вже неможливо.
Тому при одинарній точності: діапазон зсунутого порядку 1...254. Якщо виберемо зсув 128: діапазон незсунутого порядку —127...126, мінімальне нормальне число 2↑-127, максимальне дещо менше ніж 2↑126. Ой, не те. Якщо 127: діапазон незсунутого порядку —126...127, мінімальне нормальне число 2↑-126, максимальне дещо менше ніж 2↑127 і сюди входить 2↑126. Вже нормально. Аналоґічно для решти форматів. (Примітимо, що у списку з quadibloc зсув записаний як 126. Це запис того ж зсуву, але для fraction view: 0.5⩽M<1. Для left units view (1⩽M<2) те ж саме буде дорівнювати 127, як в нашому розрахунку. Здається, в сторінці з quadibloc там все підлаштовано під fraction view. Я не перевіряв всі ≈60 форматів. Але так легко заплутатись...) У стандарті ця вимога виглядає зараз так: для незсунутих порядків нормалізованих чисел, e_max + e_min = 1 (розділ 3.3). Видно, що це передбачає парну кількість можливих значень порядку.
Перевірити своє розуміння зсуву порядку дуже легко. У будь-якому з стандарних транспортних форматів, при будь-якій довжині представлення, бітова послідовність 0100...00 має інтерпретуватись як +2.0 (а 1100...00, відповідно, —2.0).
Приклад представлення для деяких чисел (обидві точності): показуємо побітно, як «sign : biased exponent : significand»:
1: single: 0:01111111:00000000000000000000000 double: 0:01111111111:0000000000000000000000000000000000000000000000000000 -1: single: 1:01111111:00000000000000000000000 double: 1:01111111111:0000000000000000000000000000000000000000000000000000 1.5: single: 0:01111111:10000000000000000000000 double: 0:01111111111:1000000000000000000000000000000000000000000000000000 2: (дивитись вище про перевірку розуміння поля порядку) single: 0:10000000:00000000000000000000000 double: 0:10000000000:0000000000000000000000000000000000000000000000000000 3.141592653589793: single: 0:10000000:10010010000111111011011 double: 0:10000000000:1001001000011111101101010100010001000010110100011000 +Infinity: single: 0:11111111:00000000000000000000000 double: 0:11111111111:0000000000000000000000000000000000000000000000000000 +qNaN (наприклад, впишемо якийсь візерунок у payload): single: 0:11111111:10001110000111100000111 double: 0:11111111111:1000111000111000011110000111100000111110000001111110
Ну, читати це очами не має сенсу. Але можна погратись, модіфікуючи число на маленький доданок і дивитись, що отримуємо.
Що ще важливо (або цікаво, як кому) у представленнях згідно з IEEE754?
1) У будь-якого числа, що взагалі передається у якомусь форматі, є тільки одне представлення (юрист додає: в цьому форматі) — це витікає з фіксації цілої частини мантиси.
2) Якщо дивитись на бітове представлення як на ціле число, то послідовність представлень від +0 до максимального скінченного числа і наступним кроком до +∞ (і таким же чином від —0 до -∞) відображається в неперервну монотонну послідовність таких цілих чисел.
Це дійсно зручно. По-перше, сортування за такими числами, якщо треба: прості нечислові маніпуляції — і з числа можна зробити ключ для сортування у стилі бази даних, чи взагалі як рядок бітів (про це буде в частині 6). По-друге, отримання сусідніх чисел для пошуку ефектів нестійкости або покращення результатів (C: функції nextafter, nexttoward).
Одинарна і подвійна точність, і не тільки. Чому дві і в чому різниця.
Чому два різних рівня точности (зі стандарта 1985) або навіть більше?
Якщо звернутись до того ж списку форматів, видно, що майже всі архітектури мають два рівня точности. Наскільки це застаріле? Насправді не дуже: витрати часу на подвійну точність десь в 5 разів більше ніж на одинарну у сучасних відеокарт NVidia (стверджують в Інтернеті), близькі цифри будуть у інших архітектур (враховуючи і час на передачу даних, і складність самих обчислень). Тобто якщо можна скоротити точність і пришвидшити обчислення, то цим треба користатись: буде більше FPSів, чи складний розрахунок виконається скоріше...
Почнемо з порядку. Вимога є такою: формат з одинарною точністю має підтримувати всі основні фізичні константи. Їх багато, але до них входять з максимальним (за модулем) порядком такі константи, як число Авоґадро (десятковий порядок +23) і константа Планка (порядок —34). Якщо ми виділяємо 7 біт на порядок, його діапазон буде приблизно ±64, а само число — приблизно (до порядку) 10↑±19; це мало. Якщо виділяємо 8 біт, діапазон порядку буде приблизно ±128, а само число — приблизно 10↑±38. Цього достатньо, наступний рівень (десятковий порядок ±76) вже є перебором. Скільки біт буде виділено на мантису (significand в стандарті) — решта місця в форматі.
Далі, мантиса: скільки цифр треба підтримувати для одинарної точности? Це вже суто емпіричні результати, але їх має бути не менше 4, краще 5 — щоб з такими числами мало б хоч якийсь сенс працювати. Насправді, IEEE854 каже: «b↑(p-1) ⩾ 100000» (p — кількість цифр вибраної основи, b — основа); для основи 10 це 6 цифр (далі буде, що для
І ось тут починається про що подумати... Скільки десяткових цифр можна представити, маючи 24 біти точности (включаючи прихований старший біт), це параметр min10d з таблиці? Формула така: floor((L-1)×log10(2)), при L=24. Без floor отримуємо 6.92... якщо було б в мантисі хоча б на 1 біт більше, було б 7.22 — тобто завжди представляло би 7 десяткових цифр. А так — тільки 6 цифр, і 8.589973e9 не можна представити: при конверсії во внутрішній формат і назад буде 8.589974e9 (це найменше з таких значень, більших 1; їх, звісно, ще є). Наскільки це важливо? Вже казав: коли у нас всі дані і вхідні і вихідні двійкові, це неважливо, якість розрахунків відповідає кількости двійкових розрядів. А ось якщо нас змушують все оцінювати в десятковій... така втрата одної цифри «на рівному місці» це дещо образливо. Але якщо б ґарантувалось 7 десяткових цифр, то діапазон порядку був би 10↑±19. Мало місця...
Ну, це не перший випадок. Вже 60 років тому назад була проблема при переході від IBM 709/7090/ітп. — 27 двійкових розрядів (у повному слові на 36 бітів) до IBM 360, де на мантису було не просто 24 біти — а ще й с основою 16 (6 таких цифр), для значення 0.1 залишався 21 значущий біт... фактично це були ті ж ґарантовані 6 цифр, але вже без того, що 7 у більшости випадків (формула дає 6.02). Де 7090 давала 7 ґарантованих знаків, 360 давала тільки 6. Перехід був, за деякими відгуками, катастрофою для розрахунків (не тільки через одну цифру, а, головним чином, і через різке підвищення нестабільности похибки). Подвійна точність частково рятувала, але ціною часу і місця на числа; тільки розвиток електроніки виправдав цю зміну. (Але біля 2000 року IBM і в лінію нащадків S/360, тобто SystemZ або просто Z, ввела реалізацію за IEEE754.)
Що для подвійної точности? Цей термін не означає, що точність (всі її параметри) якось стали в 2 рази більше. Фактично діють принципи:
1) Подвійний розмір повного представлення. Тотальна двоїчна ієрархія розмірів, яка закріпилась остаточно разом із
2) Порядок розширюється... тут вимога: формат подвійної точности має представляти
(Для наступного переходу до 128 біт тут розширення на 4 біти, а не на 3, тобто не
(IBM S/360 цю вимогу не виконувала, діапазони порядків були однакові. Ну так у них це було не єдиним дивним.)
Ну і на мантису знову решта. 53 біти (включаючи прихований) замість 24 це не в 2 рази, це трохи більше.
Зараз багато хто про одинарну точність і не знає. В ігробудуванні — знають, там навіть одинарна іноді надмірна, для ґрафіки використовують половинну (16 біт). В машинному навчанні, в нейромережах бувають навіть
Я сам не бачив реальних випадків потреби в точності більш ніж подвійна (як, наприклад, 128 біт), але чув про них. Реальна частка випадків використання такої точности настільки мізерна, що навряд чи буде більше 100 місць на всю планету. Навпаки, скорочена точність (менше «одинарної») це зараз типова вимога і для неї є апаратна підтримка в багатьох архітектурах. (Про це буде окремо.)
І все одно — навіщо концентруватись на двох варіантах точности? А тому, що оптимізують на те, що є «прості» розрахунки, для яких точности біля 6 десяткових цифр достатньо, і інші, для яких треба більше... і тут «раптово» виявилось, що точність більше, приблизно,
IEEE754 специфікує зараз будь-яку точність, у кроках по 32 біти представлення, універсальним чином. Формули там дивні, на перший погляд: якщо у числа N біт повного представлення, у нього має бути round(4*log2(N))-13 біт порядку. Але це є, виражене таким формальним чином, описане вище правило «до степеню 8 у наступній подвоєній точності» з розширенням з 8 до 16 (бо 4 біта зручніше, ніж 3).
Виконання операцій
Стандарт визначає обовʼязкові операції (арифметика, конверсії, аналіз значення, деякі супутні операції) і рекомендовані (більш серйозні математичні: степені, лоґарифми, триґонометрія). Всі, окрім явної конверсії форматів, на виході мають такі ж типи, як на вході (наприклад, стандартна подвійна точність).
На всі операції, що специфіковані у стандарті, є одна принципова вимога: результат має бути таким, якби він виконувався з нескінченними точністю (і діапазоном порядку — про це часто забувають) і вже фінально урізають (з округленням, заміною на 0 чи ∞) при розміщенні у вихідному операнді.
Пояснювальний набір прикладів на це краще надати в основі 10. Уявімо, що у нас 6 значущих цифр. Це значить, що ми можемо представити 120456 (як 1.20456e5), 0.000120456 (як 1.20456e-4), 120456000 (як 1.20456e8), 120450, 12045, 0.00012045, 120400, 12040, 1204... але не 1204567. Всі нулі на початку і на кінці запису не значущі, залишаємо тільки той ланцюжок, де перша і остання цифри не 0, і вже його довжина обмежена точністю мантиси.
123000 (1.23e5) + 456 (4.56e2) == 123456 (1.23456e5), результат точний.
0.123 (1.23e-1) + 0.000456 (4.56e-4) == 0.123456 (1.23456e-1), результат точний.
123000 (1.23e5) + 4.56 (4.56e1) == 123045.6 як з нескінченною точністю (що він також дорівнює 123045.6000...00... всі вже мають знати), але це 7 значущих цифр; треба замінити на щось близьке (далі буде формалізація), що можна представити. Це і зветься «округлення». Два кандидати — 123045 і 123046, який вибрати — буде далі (зазвичай, 123046).
1/3 = 0.(3) = 0.333.... (нескінченно цифр), маємо урізати і округлити, замінюємо на 0.333333 (3.33333e-1).
2/3 = 0.(6) = 0.666.... (нескінченно цифр), маємо урізати і округлити, замінюємо на 0.666667 (6.66667e-1); чому не 0.666666? А саме тому, що округлення за замовчуванням (див. далі) тут спрацювало «доверху», 0.666667 ближче.
1/7 = 0.(142857), зберігаємо 1.42857e-1 == 0.142857.
4/7 = 0.(571428), зберігаємо 5.71429e-1 == 0.571429.
Ну і так далі.
І ось з цею вимогою — якби нескінченна точність і потім урізання точности з округленням — починається... спочатку це був шалений ліс граблей для реалізаторів і для тих, кому треба швидко. Зараз, скоріше, такий собі гайок, бо все пройшли, стежки розмічені...
Округлення
Нащо потрібне округлення — більш-менш очевидно. Якщо у нас значення, наприклад, 222.9, а треба записати близьке ціле число, то 223 дасть меншу похибку, ніж 222. Округляти вчили в школі (ще вчать?), у десятковому варіанті (зазвичай, якщо найбільша втрачена цифра 5 або більше, то до більшого за модулем, інакше до меншого). Якщо значення 2.7182818284590452..., а треба зберегти тільки 7 цифр всього, залишаємо 2.718282 (останню підвищили), а якщо 8 цифр — 2.7182818 (остання залишилась незмінною).
За простою ідеєю стоїть ціла купа подробиць, які іноді стають критичними.
По-перше, визначення округлення. Інтуїтивно зрозуміле стає проблемним, коли ми намагаємось сформувати чітки правила. Округлення потрібне, коли ми не можемо напряму записати число в дозволеній формі. Ось ми можемо навести тільки 3 значущі цифри, а їх у нас 5. Потрібно вибирати з варіантів. Скільки цих варіантів і які вони? Так і хочеться сказати «два найближчих». Вже тут засідка. Ось вичислили 998+5=1003, а треба залишити тільки 3 цифри. Два найближчі значення, що мають до 3 значущих цифр, це 999 і 1000. А вибирати, насправді, треба між 1000 і 1010. Переформулюємо в коректний вигляд: одне найбільше з не більших за модулем (TZ — toward zero), друге найменше з не менших (AZ — away from zero). Або можна не за модулем, тоді треба буде далі постійно уточнювати деякі формулювання, а цього не хочеться. Якщо TZ = AZ, представлення точне, округлювати не треба. Інакше — треба вибрати щось одне.
Далі, які режими. Те, що в школі вчили — де дивимось на найбільшу втрачену цифру (0÷4 — усікаємо, 5÷9 — додаємо 1 до того що залишається), зветься «round to nearest, ties away» («до найближчого, хилиться від нуля»), або скорочено «half-away» (всі half-${something} це скорочення для round to nearest, ties ${something}). Воно корисне для фінансів (бо вже такі правила), але не для обчислень, що йдуть від фізики або статистики (ну хіба що для фінальних результатів розрахунків).
Десятиліття розвитку IT індустрії привели до вибору, як основного в наукових розрахунках, варіанту «round to nearest, ties to even» («half-to-even»), воно ж банкірське округлення, воно ж голандське (Dutch), воно ж округлення Гаусса, воно ж ще декілька назв... Наприклад, якщо до цілого, то 3.1, 3.4 округлюються до 3; 3.5, 3.6, 3.9, 4.1, 4.4, 4.5 до 4; 4.6, 4.9, 5.1, 5.4 до 5; а от 5.5 вже до 6... Чому саме такий режим?
По-перше, він гарніше компенсує різнонаправлені похибки, ніж багато інших, бо нема переваги округленню в один бік (як від 0 для half-away); тому воно було вигадане ще за декілька століть до перших компʼютерів. (І можна замислитись, чи не гарно було б його стандартизувати і у нас для фінансів.)
По-друге, перевага парним числам, як наслідок, то круглим числам (10 центів, а не 11). Це має ще один позитивний результат: оскільки нуль — парний, цей режим краще справляється з хвостовими нулями при імпорті і експорті.
По-третє, воно (разом з «half-to-odd») має гарну властивість: у послідовності операцій типу x+y-y+y-y... значення стабілізується і нема постійного зміщення на кожній парі таких операцій. Кому цікаво, чому це так важливо — дивіться книги, але в двух словах — полегшує оптимізацію.
Нарешті, воно з серед усіх «до найближчого» не суттєво складніше інших для апаратної реалізації.
А от якщо напряму викликати конверсію до цілого значення (як int() в Python, Java, C#, etc., (int) в C), відомі реалізації просто «відрізають» (truncate) дробову частину. Це також округлення — але вже в напрямку нуля (чи просто «до нуля», toward zero). Для C, C++, багато ще де, така поведінка конверсії встановлена їх стандартами.
Ці два режими — half-to-even (для всіх операцій) and toward zero (при чому тільки при округленні до цілого) — складають ґарантований набір у переважній більшости реалізацій плаваючої крапки; в багатьох ви і не знайдете іншого. Може, і не треба інших? (Звісно, тільки для основи 2)
Фактично, найбільш часто можна побачити 5 основних режимів.
1. Round to nearest, ties to even (округлення до найближчого, хилиться до парного), або half-to-even. Воно завмовчки в IEEE754 там, де режим округлення регулюється, і єдине там, де не регулюється (відходячи від стандарту), коли стандарт задовільняється в більшости інших деталей. Приклади були десятком абзаців вище.
2. Round to nearest, ties away (округлення до найближчого, хилиться в напрямку від нуля), або half-away. Має бути доступним за IEEE754 для основи 10. Він же прийнятий у нас у більшості випадків для фінансів: дивимось на першу відкинуту цифру (що не влізла в доступне місце в результаті), якщо вона від 0 до 4, округлюємо в напрямку нуля, інакше (від 5 до 9) — від нуля. Мабуть, зручно тому, що не треба визначати, після 5 ще йдуть якісь ненульові цифри, чи ні.
На тому ж наборі при округленні до цілого: 2.1, 2.2, 2.4, 2.499 заміняються на 2; 2.5, 2.50001, 2.51, 2.6, 2.9,..., 3.4, 3.49, 3.49999 заміняються на 3; 3.5, 3.50001,... заміняються на 4, і так далі. Тобто різниця порівняно з half-to-even буде при значеннях 0.5, 2.5, 4.5...
3. Round toward zero (округлення в напрямку нуля), або просто усікання дробової частини. Підтримка вимагається в IEEE754. Ще trunc() функція POSIX. Крім того, встановлено для виконання при явній конверсії в ціле, щонайменше, в C, C++, Java, C#, Python.
4. Round toward minus infinity, toward -∞ (округлення в напрямку мінус нескінченности). Підтримка вимагається в IEEE754. Ще floor() функція POSIX.
5. Round toward plus infinity, toward +∞ (округлення в напрямку плюс нескінченности). Підтримка вимагається в IEEE754. Ще ceil() функція POSIX.
Інші теоретично можливі без знущання над совою і ґлобусом: away from zero, half-to-odd, half-floor, half-ceil... — вимагаються у настільки специфічних ситуаціях, що їх не стандартизують. Хоча см. далі.
Режими 3÷5 звуться направленими (direct), тобто не до найближчого. Вони необхідні для підтримки інтервальної арифметики — якою перевіряють можливі діапазони відхилень для розрахунків: виставляється режим, проводиться така ж порція обчислень, як і в звичайному режимі, і оцінюється відхилення результату (це дуже примитивний опис, але розписувати детальніше вимагатиме специфічної математики). Це нечаста задача. Тому ви в стандартних засобах дуже багатьох мов і їх рантаймів — як Python, Java, C#, JavaScript... — не знайдете штатного методу змінити режим округлення для основних операцій. Якщо є (java.math.RoundingMode), то для специфічного використання (BigDecimal).
Які ще бувають варіанти округлення? Ось, наприклад, НБУ обіцяло з осені 2025 видаляти монети у 10 копійок. Після цього найменша грошова одиниця — 50 копійок (чи шаґів, якщо перейменують). Ціни при розрахунку готівкою мають бути округлені до 50 копійок... цей режим не відповідає жодному з описаних. Як це буде зроблено? Простіше за все домножити на 2, округлити до гривень (half-away, бо фінанси) і поділити назад на 2. Але є і інші методи. Ну, подивимось, як воно буде специфіковано. Нам тут важливе те, що можуть бути дуже неочікувані варіанти... але в залізі їх не чекайте, треба емулювати проґрамно.
Але це все не включає один режим, який рідко згадують, але він лежить в основі всіх алґоритмів обчислень: так званий rounding to prepare for shortest precision (позначимо як RPSP), він же — округлення фон Неймана (цей розумний дядько винайшов далеко не тільки базову архітектуру компʼютера;)) Що він такеє? Є проблема, яка зветься подвійним округленням (double rounding); пояснимо її простим десятковим прикладом. Ось є число 5.49 (вважатимемо, основа 10, представлення точне). Режим округлення — half-to-even або half-away, цей приклад спеціально підібраний, щоб було одинаково. Якщо ми одразу округлимо до цілого, отримаємо 5, бо 5.49 ближче до 5, ніж до 6. Але, якщо будемо спочатку округляти до десятих (0.1), а потім до цілих, отримаємо 5.5, яке для обох half-to-even і half-away округляється до 6. Ууупс, неочікувано... В статті Rounding в Вікіпедії надано приклад судового процесу в США, в якому фінальний вирок був — позбавитись ефектів такого подвійного округлення і рахувати так, щоб було еквівалентно тому, якби одразу округлювалось до фінальної точности.
І ось тут на поміч приходить RPSP. Його принцип: якщо значення хвосту, що відрізається при округленні, інші, ніж 0 і ½ від ваги найменшої цифри, що залишається, то ці 0 і ½ не можуть бути вибрані, відповідний кандидат з двох (TZ і AZ) забороняється. Якщо ми 5.49 округлюємо до десятих, то між варіантами TZ=5.4 і AZ=5.5 ми повинні вибрати 5.4, незважаючи, що 5.5 був найближчим: 5.5 дозволено тільки коли б число до того було рівно 5.5. Аналоґічно, 5.5000001 має бути округлено до 5.6 а не 5.5; 5.0001 до 5,1, а не 5.0... Головне: після цього наступний етап (як округлення до цілого в даному прикладі) буде виконано коректно незалежно від режиму цього округлення — для всіх тих, що хоч раз згадувались як «природні» в цьому циклі статей, і не тільки: toward zero, away from zero, toward -∞, toward +∞, half-to-even, half-to-odd, half-floor, half-ceil, half-to-zero, half-away, і сам RPSP. Можна розширити це до трьох кроків округлення і так далі, поки витримується, що останнє з них — за цільовим режимом, а всі попередні — в RPSP, тобто сам RPSP теж можна додати до цього списку.
Виконання того ж методу для основ 2, 8, 16 — див. нижче.
Умова коректности застосування RPSP — між кожними кроками округлення повинно даватись, що попередній результат зводиться до одного 4 варіантів (хвіст точно 0, менше ½, точно ½, більше ½ від ваги цифри перед ним), що для основ 8, 10, 16 не накладає допоміжних умов (можна на одну цифру), а для основи 2 — щоб був крок точности між попереднім і наступним округленням не менше ніж 2 біти (і тому не менше 2 біт «хвосту» перед фінальним округленням).
Перша згадка про цей метод, наскільки я накопав, із 1949 року, з робіт над EDSAC (для двійкових чисел). Але його рідко згадують саме з іменем фон Неймана, і це мене дивує. (Цікаво, що в подальшому сам Джон фон Нейман виступав проти рухомої крапки, вважаючи, що фіксованої достатньо тому, хто правильно спроєктував проґраму. Але не забуваємо, що це ще початок
І ось цей метод фактично використовується в усіх (що я бачив... а нахіба робити інакше?) сучасних реалізаціях плаваючої крапки: заданий контекстом режим округлення використовується тільки на фінальній стадії. Основні тіла виконання операцій, як додавання, множення і т.п., вираховують проміжний результат до тої точности, що достатня для виробки проміжного результату, який отриманий при безальтернативному RPSP. Маємо перше округлення — при усіченні хвоста (потенціально, нескінченного) і друге — при формуванні кінцевого результату.
Наприклад: нехай у нас мантиса 6 десяткових цифр. Складаємо 123000 з 4.56. Проміжний результат на 7 цифр (бо треба ще одну як раз) з точного 123004.56 скорочується до 123004.6. Далі працює фінальне округлення: half-що-завгодно і toward +∞ дадуть 123005; toward zero і toward -∞ дадуть 123004.
В десятковому варіанті 4 описані варіанти хвоста мапляться: 0; 1÷4; 5; 6÷9 тощо. Коли відрізається зайва частина хвоста значення, цифри крім 0 і 5 зазвичай не модифікуються, бо не треба; але 0 і 5, якщо відрізана частина була не равна 0, заміняються на 1 і 6 відповідно.
В двійковому варіанті треба 2 біти на хвіст. І вони натурально «мапляться» на операцію, що зветься jamming shift: якщо в тому, що не попало навіть у хвіст, є ненульові біти, наймолодший біт збереженого хвоста ставиться в 1 (могло бути, що там вже 1, тоді він не змінюється). Кому цікаво в максимальних подробицях — поґуґліть рахувалочку «guard, round, sticky». Те ж саме можна робити при основах 8 і 16, коректуючи тільки наймолодший біт. Зветься воно в такому випадку округленням до непарного, to odd (не плутати з округленням до найближчого і схилом до непарного, half-odd! це різні).
Іноді в бібліотеках RPSP є в явному вигляді: якщо бачите назву режиму типу ROUND_05UP, це воно. Дивує, що вважають за доцільне дати йому назву, специфічну для десяткової арифметики.
(Найповніший набір в залізі, що я десь бачив, це в IBM SystemZ. Направлені (direct): до нуля (toward zero), від нуля (away from zero), до найближчого меншого (toward -∞), до найближчого більшого (toward +∞). Напівнаправлені (вибір тільки посередині): half-to-even, half-away (від нуля), half-toward-zero. Нарешті, RPSP. Навіть з їх специфікою я не впевнений, що це все колись комусь реально потрібне в залізі.)
Але дивно, що RPSP майже всюди не дають за межі внутрішньої кухні, хоча фактично він є всюди.
Стандарт IEEE754-2008 вже вимагає присутности окремої функції перетворення до цілого для кожного з режимів округлення, і загального режиму в контексті (точніше, двох — роздільно для двійкових і десяткових чисел — якщо десяткові підтримуються) для решти операцій. Але з нього уже почалась рекомендація давати можливість задавати режим не контекстом, а в даних самої операції. Поки що я не бачу в засобах популярних мов навіть явного режиму для перетворень. Краще, що є — виставлення режиму в контексті і його застосування. Ось для C/C++ маємо:
- Поточний режим:
rint,lrint*,llrint*,nearbyint*. (Всюди тут де «*» - без суфіксу для double, з `f` для float, з `l` для long double; що таке long double — окрема тема.) - До нуля: пряма конверсія типу, trunc*.
- До -∞:
floor*. - До +∞:
ceil*. - Half-away:
round*. - Загальний контекст для інших операцій (не до цілого) і там де поточний режим:
fesetround()і `#pragma STDC FENV_ACCESS`, проте не всі компілятори і рантайми їх чесно підтримують.
Якось нерівно. Ну, стандарти системи «лебідь, рак і щука» вони завжди такі.
У розробників процесорів теж інерція, цілком зрозуміла: проблеми сумісности і леґасі у них тяжче, ніж в софті. В RISC-V додали можливість задавати режим округлення в кожній команді (свій, якщо не посилається на дінамічний, що з реґістру контексту). Це єдина серйозна ISA (Instruction Set Architecture), що стартувала з нуля після ≈2005 року. Ну, поки у неї активний розвиток, подивимось на перспективи.
Що тут можна сказати як висновок?
1. Спостерігаємо тяжкий розрив між типовими вимогами практики і рекомендаціями загального стандарту, який розрахований на значно більш математичні задачі, ніж в переважній практиці.
2. Треба завжди тримати на увазі, яке округлення нам потрібне, і уважно слідкувати за деталями. Ну як завжди.
Виконання простих операцій: базовий алґоритм
Один з найвідоміших в opensource світі и найдоступніший зараз варіант проґрамної реалізації базових операцій (арифметика) з числами з рухомою крапкою це Berkeley SoftFloat; воно підтримується зараз командою, близькою до розробки RISC-V. На цій бібліотеці краще всього дивитись на принципи роботи, звісно, роблячи поправки на теорію.
(Щонайменше, є ще soft-fp в glibc і її копія в libgcc. Але там складно щось роздивитись через макри, що на макрах сидять і макрами погоняють.)
Візьмемо додавання і віднімання
При додаванні, менше за порядком число зсувається праворуч, але цей зсув особливий: він зветься jamming shift (повторюю, що було декілька абзаців тому), при ньому має бути визначено, чи у втраченій частині був хоч один біт не 0; якщо був, найменший в тому, що залишилось, ставиться в 1. Якщо подумати, легко видно, що це варіант RPSP (см. частина 1) для основи 2 — те саме округлення до непарного. Дані з акумуляторів складаються, і хвіст аналізується для вибраного фінального округлення на 4 випадки (точно 0, менше ½ від одиниці мінімальної цифри до хвоста (voting digit), точно ½, більше ½); можлива ще корекція порядку, якщо округлення привело до переносу в новий найбільший розряд; нарешті, результат пакується в вибраний формат результату.
І знову — краще показувати на десяткових розрахунках, хоч у них і дещо специфіка. В цих прикладах буде 6 десяткових цифр фінальної точности (і 8 проміжної, див. деталі), обовʼязкова нормалізація і незсунутий порядок в right units view, всі числа додатні, при відніманні перше більше другого (інакше треба слідкувати за знаками).
Операція додавання: в усіх акумуляторах 1 десяткова цифра перед (на перенос в наступний розряд) і 1 після основного значення (на округлення), відділяємо їх знаком ’|’. Те, що втрачаєтся при зсуві, відділятиметься знаком ’‖’. Загальний алґоритм:
— розібрати значення на порядок і мантису, розмістити вхідні дані в акумуляторах;
— зрівняти порядки: спочатку, якщо треба, зсув того, що з більшим порядком, вліво (даючи більше місця для точности; максимум зсуву — зрівняти порядки з іншим доданком), потім зсув вправо того числа, що з меншим порядком (може вимагати проміжне округлення за RPSP);
— скласти значення в акумуляторах;
— перевірити на перенос у наступний разряд (може викликати ще один зсув вправо на 1 цифру з проміжним округленням);
— фінальне розміщення результата (і тут впливає вибраний режим округлення).
Доданок 1: 123 == 123000e-3; акумулятор 1: 0|123000|0;
Доданок 2: 456 == 456000e-3; акумулятор 2: 0|456000|0;
Порядки однакові, складаємо: 0|579000|0; округлення немає; записуємо фінальне число 579000e-3, тобто 579.
Доданок 1: 123000 == 123000e0; акумулятор 1: 0|123000|0;
Доданок 2: 456 == 456000e-3; акумулятор 2: 0|456000|0;
Доданок 2 має бути зсунутий вправо на 3 цифри: 0|000456|0; тепер це фактично представляє 456e0, а не 456000e-3, хоча значення не змінилось;
Порядки тепер однакові, складаємо: 0|123456|0; округлення немає; записуємо фінальне число 123456e0, тобто 123456.
Доданок 1: 123000 == 001230e2; акумулятор 1: 0|001230|0;
Доданок 2: 456 == 456000e-3; акумулятор 2: 0|456000|0;
Доданок 1 має бути зсунутий вліво на 2 цифри: 0|123000|0; тепер це фактично представляє 123000e0, а не 001230e2;
Доданок 2 має бути зсунутий вправо на 3 цифри: 0|000456|0; тепер це фактично представляє 456e0, а не 456000e-3;
Порядки тепер однакові, складаємо: 0|123456|0; округлення немає; записуємо фінальне число 123456e0, тобто 123456.
Доданок 1: 123000 == 123000e0; акумулятор 1: 0|123000|0;
Доданок 2: 45.6 == 456000e-4; акумулятор 2: 0|456000|0;
Доданок 2 має бути зсунутий вправо на 4 цифри: 0|000045|6;
Порядки тепер однакові, складаємо: 0|123045|6; треба робити округлення. Дивимось на режим (з 5 вказаних в стандарті):
Якщо half-to-even, half-away, toward-plus-infinity: додаємо 1 до основної частини, проміжний результат 0|123046|0; фінально 123046e0;
Якщо toward-zero, toward-minus-infinity: не додаємо 1, проміжний результат 0|123045|0, фінально 123045e0.
Доданок 1: 123000 == 123000e0; акумулятор 1: 0|123000|0;
Доданок 2: 4.56 == 456000e-5; акумулятор 2: 0|456000|0;
Доданок 2 має бути зсунутий вправо на 5 цифр: 0|000004|5‖6, хвіст, що втрачається (після «‖») ненульовий, тому згідно з RPSP міняємо найменш вагому з залишених — 5 на 6, результат 0|000004|6;
Порядки тепер однакові, складаємо: 0|123004|6. Округлення — дивимось на режим:
Якщо half-to-even, half-away, toward-plus-infinity: додаємо 1 до основної частини, проміжний результат 0|123046|0; фінально 123005e0 == 123005;
Якщо toward-zero, toward-minus-infinity: не додаємо 1, проміжний результат 0|123004|0, фінально 123004e0 == 123004.
Доданок 1: 789012 == 789012e0; акумулятор 1: 0|789012|0;
Доданок 2: 654321 == 654321e0; акумулятор 2: 0|654321|0;
Складаємо: 1|443333|0; маємо зсунути ще раз вправо на 1: 0|144333|3; округлення за RPSP, не треба модифікувати найменшу цифру, що залишилась.
Якщо half-to-even, half-away, toward-zero, toward-minus-infinity: не додаємо 1, проміжний результат 0|144333|0, фінально 144333e1 == 1443330.
Якщо toward-plus-infinity: додаємо 1, проміжний результат 0|144334|3, фінальний 144334e0 == 1443340.
Фінально найскладніший для додавання (і перенос, і коректування після зсуву):
Доданок 1: 999987 == 999987e0; акумулятор 1: 0|999987|0;
Доданок 2: 68.4321 == 684321e-4; акумулятор 2: 0|684321|0;
Необмежено точний результат дорівнює 1000055.4321;
Доданок 2 має бути зсунутий вправо на 4 цифри: 0|000068|4‖321; найменш вагома залишена дорівнює 4, тому не модифікується незалежно від втраченого хвоста;
Складаємо: 1|000055|4; нам треба ще один зсув вправо на 1 за RPSP, отримуємо спочатку 0|100005|5‖4, при зсуві втратили не 0, тому 5 міняємо на 6; фінально зараз 0|100005|6;
Якщо half-to-even, half-away, toward-plus-infinity: додаємо 1 до основної частини, проміжний результат 0|100006|0, фінальний 100006e1 == 1000060;
Якщо toward-zero, toward-minus-infinity: не додаємо 1, проміжний результат 0|100005|0, фінальний 100005e1 == 1000050. (Вірно, оскільки 6 цифр мантиси, то крок значень тут 10.)
Інші арифметичні операції робляться дуже подібно, але зі своєю специфікою. При відніманні (і множенні) можливо, що порядок зменшиться на 1 при вимозі тримати повний хвіст на округлення, тому треба на 1 біт більше в хвості — щонайменше 3; цей біт зветься guard bit. (Для десяткової це не 3 біти, а 2 цифри.) Чому 3 біти або 2 десяткових цифри? Візьмемо випадок: 1000000 — 3.61.
Помилковий контрприклад: якщо рахуємо з 1 цифрою і рештою алґоритму, як для додавання:
Зменшуване: 1000000 = 1.00000e6; акумулятор 1: 0|100000|0;
Відʼємник: 3.61 = 3.61e0; акумулятор 2: 0|361000|0;
Точний результат без округлення — 999996.39;
Відʼємник має бути зсунутий вправо на 6 цифр: 0|000000|3‖61; залишається 0|000000|3;
Різниця 0|099999|7; коректуємо зсувом, 0|999997|0 (зсунули вліво на 1 цифру, бо найстарша була нуль); відповідь 999997? Вона невірна, при half-to-even, half-away має бути 999996.
Чому втратили? Тому, що різниця порядків зменшуваного і відʼємника більше 1, і ще виник ведучий нуль.
Змінимо на наявність двох цифр хвоста після основного значення (а резервна в голові не потрібна, ми віднімаємо). Повний алгоритм:
Операція віднімання: в усіх акумуляторах 2 десяткових цифри після після основного значення (на округлення), відділяємо їх знаком ’|’. Загальний алгоритм:
— розібрати значення на порядок і мантису, розмістити вхідні дані в акумуляторах;
— зрівняти порядки шляхом зсуву вправо того числа, що з меншим порядком (може вимагати проміжне округлення за RPSP);
— відняти значення в акумуляторах;
— перевірити на втрату порядку (може викликати зсув вліво);
— округлення (може давати перенос назад у старшу цифру і тому зсув вправо на одну цифру, але при цьому найменша дорівнює 0, тому можна не виділяти на неї місце);
— фінальне розміщення результата (і тут впливає вибраний режим округлення).
Зменшуване: 1000000 = 1.00000e6; акумулятор 1: 100000|00;
Відʼємник: 3.61 = 3.61e0; акумулятор 2: 361000|00;
Відʼємник має бути зсунутий вправо на 6 цифр: 000000|36‖1, залишається 000000|36;
Різниця 099999|64; тепер нормалізуємо зсувом на 1 вліво, отримуємо 999996|40. Тут нарешті застосовуємо вибране округлення, і маємо коректний результат після округлень:
Якщо half-to-even, half-away, toward-zero, toward-minus-infinity: 999996e0 == 999996.
Якщо toward-plus-infinity: 999997e0 == 999997.
(Коректність розширення лише на 1 цифру вимагає доведення, але воно тривіальне. Читачі можуть зробити це самі.)
І ось тут можна прочитати головне, що інтернет-пошук видає по словах «guard round sticky», де те ж саме описується, по більшости, для двійкової арифметики. На основі 2 принципи ті ж, звісно, але маніпуляції з хвостом виглядають інакше. При такому ж запису, якщо 2 біти (додавання і ділення),
— xxx|00‖Z це |00‖ перетворюється на 01, якщо в Z (скільки завгодно бітів) хоча б одна одиниця, інакше залишається 00;
— xxx|01‖Z це |01‖ завжди залишається 01;
— xxx|10‖Z це |10‖ перетворюється на 11, якщо в Z хоча б одна одиниця, інакше залишається 00;
— xxx|11‖Z це |11‖ завжди залишається 11;
те ж саме я описував раніше як 4 варіанти 0, 0<y<½, ½, ½<y<1. Постійний привіт від дядьки Яноша (Джона).
З аспекту роботи над найближчими цифрами хвоста, множення принципово не відрізняється від віднімання, а ділення — від додавання. Для множення і ділення варто зупинитись, коли отримана достатня кількість бітів частки і ми знаємо, проміжний залишок 0 чи ні (тоді пишемо ще одну одиницю, для jamming shift цього досить).
Виглядає все прямолінійно і навіть просто (да, не для того, хто перший раз читає, але якщо є мінімальний досвід). Але за цею простотою стоять таки десятиліття розробки і пошуку точного і оптимального підходу. Навіть зараз мало які джерела описують в явному вигляді саме двохкрокове округлення при арифметичних операціях і проміжний режим, проходячи це якось побіжно, а це центральне в алґоритмах. Теорія коректного округлення, звісно, була раніше ніж IEEE754 (якщо вважаємо з винаходу округлення фон Неймана, то з
Помилки і винятки
Що буде, якщо попросити стандартну бібліотеку вирахувати квадратний корінь з мінус одиниці? Ні, воно не видасть «i», бо відповідь має лежати в множині дійсних чисел. Але такої відповіді нема, і маємо щось зробити з цею помилкою. Або ми множимо максимальне число, що можна представити, на само себе; звісно, представити це вже не можна. Що робити?
Можна перелічити весь спектр варіантів, як у різних стилях і мовах обробляються помилки всіх типів задач, не тільки рухомої крапки:
1. Видається результат такого типу даних, що приймає значення або `Result(результат)`, або `Error(помилка)` (у дужках — параметр варіанта), наприклад, `Result(1.4142135)` або `Error(Domain)`. Поширений у інтерфейсах ядра (Unix, Windows обидва; помилка може доставлятись побічним каналом, як псевдо-змінна errno), у мовах, як Rust (тип результату явно є типом-варіантом). Результат треба перевірити, по компонентах або через відповідність шаблону (pattern matching), або «розгорнути» (unwrap; у випадку помилки це інша помилка і обробляється окремо).
2. Виконання завершується в цьому місті, переходимо на якийсь обробник. Багато де це зроблено через «обробку винятків» («exception handling»), включаючи її «структурований» варіант (коли є якась ієрархія винятків).
3. Ставиться якийсь ґлобальний (або свій на кожну нитку) «флаґ», «прапорець» («flag»), який треба перевіряти в ключових місцях... а що повернути з функції як результат замість числа-значення?
Можна довго казати (або не казати і сперечатись з цим), що помилки мають оброблені якомога раніше. Але не можна ґарантувати універсальну наявність обробки з миттєвим припиненням виконання. Тому стандарт зроблений під «розширений» пункт 3: і помилка детектується і підіймається флаґ, і повертається спеціальне значення «не число» (тобто, «тут нема ніякого числа») або «нескіченність» (тобто «більше за модулем, що ми могли представити, і подробиць більше немає»). А ось якщо ще є і підтримка миттєвої ґенерації винятка, то і виняток може бути запущений в обробку (і на це є окремі правила).
Є контекст рухомої крапки, який зазвичай є частиною стану процесора. В ньому є, щонайменше, режим округлення і накопичувальні флаґи помилок. Оскільки будь-який сучасний рантайм багатозадачний, контекст звʼязаний з задачею (може бути названа «процес» (process) або «тред», «нитка» (thread), або корутина (coroutine)...) і при зміні задачі, що виконується, контекст задачі, яка зупинена, має бути збережений, а тої, що починає виконуватись — піднятий зі збереженого і активований. (Тут, як завжди, свої проблеми: наприклад, чи виконує фреймворк асинхронного виконання цю вимогу для корутин? Але це за межами нашого циклу.) Існує пʼять стандартних помилок, які має підтримувати реалізація, що відповідає стандарту:
1. Переповнення (overflow): фінальний результат операції, який вираховувався як скінченний, отримав порядок, який завеликий і його не можна представити в вибраній точності. Звернімо увагу, якщо на вході вже була Infinity, ця помилка не визначається. Зґенероване значення — Infinity з потрібним знаком.
2. «Недоповнення» (underflow; ну нема кращого перекладу, беремо жарґон, він хоч легко зрозумілий). Виникає у випадку, коли результат субнормальний — менше за модулем ніж мінімальне нормалізоване число (це включає і випадок фінального нуля) і був округлений. Звернімо увагу: якщо округлення не було, помилка не ґенерується. Можна описати цю помилку як «серйозна втрата точности, за межами звичайного округлення». Зґенероване значення — найближче до потрібного, з урахуванням округлення; може бути нуль, тоді знак цього нуля показує, з якого боку підійшли.
3. Невалідна операція (invalid operation): результату просто нема для таких арґументів. Випадки, як квадратний корінь чи лоґарифм відʼємного значення, арксінус значення за модулем більше 1... Тут все очевидно. Але ще вона визначається при арґументі signaling NaN при обчислювальних операціях, при будь-якому NaN при сиґналізуючих порівняннях (буде далі), при конверсії в ціле з виходом за межі представлення... Повний список можна знайти в стандарті, він важливий для «стандарт-юристів» і реалізаторів. Зґенерованим значенням, якщо не миттєвий виняток, буде qNaN (зазвичай зараз, з канонічним payload).
4. Ділення на нуль (division by zero). Назва цеї помилки не дуже коректна; в деяких джерелах, наприклад, POSIX.1, man math_error(7) в Linux, cppreference.com, названа «pole error» («помилка полюса»?), це ближче до сути. У Microsoft було навіть «singularity error». Проте IEEE754 зберігає назву про ділення. За буквою, це випадок, коли результат — «точна» нескінченність (а не велике число, «округлене» до нескінченности, як для Overflow) при скінченних арґументах. Є тільки два випадки, коли ґенерується в операціях, визначених стандартом — ділення не-нуля на нуль і лоґарифм від нуля. В ідеалі, могло б бути також при танґенсі прямого кута, але оскільки π/2 не можна точно представити при обмеженій точності, то практично цього варіанта не існує. (За межами стандарту, в libc, додаються lgamma(), tgamma() при відʼємному цілому арґументі... тут більше варіантів.) Зґенероване значення — Infinity з потрібним знаком.
Мені не дуже зрозуміло, і не зміг знайти, навіщо відділяти таку помилку від переповнення чи невалідної операції. Як припущення, це відокремлення результатів такої операції (яка дає нескінченність при всіх режимах округлення) від, наприклад, 1e308*2 в double, яка при режимі округлення до нуля чи -∞ дасть максимальне можливе double (1.79...e308).
5. Неточний результат (inexact). Ґенерується при будь-якому округленні, коли проміжний теоретичний «нескінченно точний» результат відрізняється від того, що було представлено як результат операції при обмеженні на точність і порядок. Це включає в себе також випадок, коли скінченне значення замінюється на нескінченність. Зґенероване значення — результат округлення.
Може бути впізнано декілька помилок на одну операцію (так, Underflow завжди супроводжується Inexactʼом).
З помилкою inexact цікава ситуація: для типових математичних задач, що походять з фізики, ця помилка буде ґенеруватись майже завжди. Банальні 1/3, 1/10 вже ґенерують її. Складання таких чисел. Майже будь-які операції з числами, отриманими з реальних джерел, як сенсори. Тому там ця помилка майже не має значення. А ось за межами такої математики починається вже територія її корисности. Наприклад, якщо ви вимушені підтримувати «економічну» арифметику через рухому крапку (другий домен із початку частини 1), перевірка на inexact допоможе знайти випадки втрати точности.
Що можна робити з флаґами помилок? Там, де з ними взагалі можна працювати, підтримуються такі операції: підняти чи опустити якусь підмножину флаґів; прочитати поточний стан. (Для C: feraiseexcept(), feclearexcept(), fetestexcept().) Цього достатньо і для того, щоб детектувати помилки з попереднього коду, і щоб реґулювати, якщо треба виробити результат своєї функції, що б там в її середині ні траплялось.
І тут очікую(?) деяке виття від тих, хто проґрамує на Java, C#, JavaScript, Python і майже всьому, що виходить за рамки компільованих мов з ручним керуванням памʼяттю, що для них ці засоби (і режим округлення, і помилки-винятки) недоступні... і для чого взагалі я це згадував, коли розробка на них займає 90+% всіх витрат ресурсів проґрамістів? Не інакше, я просто дражнюсь... Хоча ні, якщо за 30+ років (Java) нікому там це не стало потрібним, то навряд чи буде і далі. Кому потрібне, викликатиме FFI переходники в C, або бібліотеки, де той же закат сонця зроблений вручну (зі втратою виробности), як decimal в Python.
Якщо вам доступний контроль цих помилок: треба чистити їх перед якимось блоком обчислень і перевіряти після цього. Звичайно це не дуже великий блок. Іноді це виконується навіть з одною операцією. Наприклад, якщо ви присвоїте значення типу float змінній типу int, а значення не можна представити в int навіть з округленням до цілого в напрямку нуля, то вихідний int міститиме щось невідоме, залежне від платформи. Тільки перевіривши контекст задачі на помилку «invalid operation», можна твердо визначити, що була проблема. (Звісно, ще є шляхи, що дорожче, як же без цього. Викликати rintf(), цілий результат перетворити назад у float, підкоректувати на 1 за вимогою, порівняти на рівність. Як же без цього. Але ціна, ціна...)
Sydorenko Oleg
Roman Gomolko
2 коментарі
Додати коментар Підписатись на коментаріВідписатись від коментарів