Плаваюча (рухома) крапка. Частини 1-2: ввідна, перші граблі

Каталог частин (поповнюваний по мірі додавання):
— Плаваюча (рухома) крапка, частина 3: Сучасний стандарт і все навколо — основи.
— Плаваюча (рухома) крапка, частина 4: Сучасний стандарт — тонкі деталі і проблеми.
— Плаваюча (рухома) крапка, частина 5: Текстовий експорт-імпорт.
— Плаваюча (рухома) крапка, частина 6: Порівняння.
— Плаваюча (рухома) крапка, частина 7: десяткова рухома крапка.
— Плаваюча (рухома) крапка, частина 8: різне. Субнормальні. Розширений 80-бітний формат. Баґ 323.
— Плаваюча (рухома) крапка, частина 9: RDBMS і SQL: земля невизначеності.
— Плаваюча (рухома) крапка, частина 10. Двійковий експорт-імпорт. Зовнішні машинні транспорти.

Навчальні нотатки (в ролях по обидва боки від кафедри). Що українською, що заліською я не бачив компактного зведення основ в потрібному обʼємі (є дуже багато, але уривками), тому вирішив звести у звʼязний текст і в процесі оновити власне розуміння. Текст дуже щільний, навіть хоча математика обґрунтувань виключена (інакше було б разів в 5-10 товстіше). Якщо критично не вистачає деталей, креслень і тп. — коментуйте. Взагалі тут, може, 1/10 від того, що можна було б розповісти, тема реально бездонна, але і того, що тут буде, більше, ніж реально потрібно більшости присутніх:)

Для того, кто більш-менш в темі, більшість змісту кожної частини це «ази». Але не пропускайте зовсім: деякі деталі можуть бути несподіваними.

У деякому сенсі, присвʼячується 40-річчю стандарту IEEE754-1985. Це навіть лякає, скільки років вже IT...

Про саму назву

Кома («1,23») чи крапка («1.23»)? Традиції різних країн дають різне. Континентальна частина Європи надає перевагу комі. Анґлія і США надають перевагу крапці. ISO дозволяє і одне, і друге. Але саме в IT фактично крапка зараз в багато раз «сильніша» за кому, бо традиції IT склались по більшости в США. Це проявляється багато в чому, від назв математичних символів (≤ або ⩽? в Unicode саме перший позначений як LESS-THAN OR EQUAL TO, хоча по Європі другий більше прийнятий), лапок (нумо доведіть якомусь Google Docs, що треба писати «так» або „так“, а не “так”«, і щоб він автоматично підміняв лапки як нам треба), і до формату запису дати. Зараз зупинимось тут, продовжимо при обговоренні текстових представлень.

В українській вікіпедії є стаття «Число з рухомою комою». Угу, там «кома», бо стандарти за межами IT у нас вимагають коми.

Чому «рухома»? Мабуть, автору терміну не сподобалось «плаваюча» (floating). В сучасній українській не дуже цінують як дієприкметники активного способу, так і прикметники, що вийшли з них (хоча «гарячий», «діючий», і т.д. — закріпились). В анґлійській чітко floating (floating-point); в російській «плавающая точка/запятая»; не буду тут сперечатись, просто маємо на увазі. В цьому циклу вони будуть використані як повні синоніми.

Головна ціль такого типу даних

Є два принципово різних домени, де використовується плаваюча (рухома) крапка, як представлення чисел. Перший — можна назвати природним, бо він 1) використовуєтся переважно для задач математичної фізики (тобто фізичних розрахунків у числах), 2) був тим, для кого цей тип представлення був створений з самого початку. В ньому переважно використовується основа 2 (см. далі). Другий — для фінансів і схожих задач; в ньому підхід з рухомою крапкою це вторинне, що походить від незадовільности наявних засобів проґрамування. В ньому переважно використовується основа 10.

Отже, почнемо для першого домена.

Ми знаємо цілі числа. Ми думаємо цілими числами у більшости випадків. Але деякі речі неможливо або, що імовірніше, немає сенсу представляти цілими числами.

Тут і далі, знак «↑» позначує піднесення у степінь; так, 2↑5 = 32, e↑(i×π) = -1.

Ось число Авоґадро — скільки молекул в одиниці речовини, в якій стільки ґрамів, скільки молекулярна маса речовини. Воно приблизно дорівнює 6.022×10²³. Ми можемо його записати в вигляді 602214082110000000000000000 (найбільш точне значення як на сьогодні), але заплутаємось у нулях в кінці — більш того, вони невірні (а вірних значень ми і не знаємо, точність не дозволяє, і навіть ті цифри, що я вказав, розбігаються в різних вимірах). Краще зрізати зайве і записати просто скорочений факт «а потім в приблизному значенні ще 16 нулів». Ну а потім перемістити крапку і сказати «23 нулі або вказані цифри» (і знову, з якоюсь точністю).

А ось константа ґравітаційної взаємодії, яка дорівнює, знову приблизно, 6.6743×10⁻¹¹ (ньютонів помножити двічі на метр і двічі розділити на кілоґрам). І знову, можемо записати у вигляді: 0.000000000066743. А нащо нам знову такі довгі ланцюжки нулів, які треба кожного разу рахувати і не збитись? Пишемо степінь числа 10 і в зручній формі. Що тут спільного, крім факту дивних чисел з багатьма нулями? Насправді є одне: фізика. Використання таких чисел росте з фізики і як її проєцують на математику.

Звісно, за межами фізики теж багато прикладів. Ми кажемо «в державних сховищах золота на 20 мільярдів доларів» чи «внутрішній борг досяг 2 тріліони беріллієвих туґриків». Практично це є те ж саме, числа узагальнюються до великих порцій, і нам вже майже неважливо, це 20 мільярдів 103 міліони 217 тисяч чи 20 мільярдів 107 міліонів 800 тисяч. Але...

Що головне в тому, де ми використовуємо числа з рухомою крапкою? Нам важлива відносна похибка значення. Ми можемо взагалі його не виміряти з точністю до одиниці (як з числом Авоґадро; можуть тут впливати квантові ефекти?), або одиниця це забагато (як з константою ґравітаційної взаємодії). Абсолютна похибка нас цікавить тут рівно настільки, наскільки вона впливає на відносну чи залежить від неї.

А ще важливіше сказати, що у типових цілях використання рухомої крапки нас в принципі влаштовує, що є якась похибка, що ми працюємо з неточними значеннями. (Тут треба зупинитись і викарбувати це твердження на камені.) Це — фундаментальна ознака, яка задає, коли, як і наскільки ми можемо використовувати такі числа.

А тепер до другого домену. В ньому все навпаки: ми не маємо права на похибку там, де не дозволили її явно.

Перш за все, це у фінансах. Бо всі правила вимагають точного представлення значень, до мінімальної одиниці (цент, пенс, копійка, чи шаґ, якщо буде введений) і ідеального сходження рівнянь. Всі суми повинні зійтись. Скільки знято з якихось бухґалтерських рахунків, стільки ж має додатись до інших. Тотальний нуль в сумі всіх, і не дай вам бог помилитись... хоча бог простить, а податкова — ні.

Наприклад, ви не можете розділити 100 гривень на трьох так щоб всім дісталось по 33.33. Ще копійку треба буде кудись діти. Або один отримає 33.34, або вона піде у комісію операції... І 200.00 не можна розділити як кожному 66.66 чи 66.67. Треба 66.67, 66.67, 66.66. Звісно, такі операції нечасті, частіше ділять на дві частини, з яких одна — комісія, або комісія береться взагалі як ще один доданок. Але і тут є проблеми — буде розказано далі.

Тому вже на дуже ранніх етапах розвитку компʼютерів світ фінансів відірвався від світу тих обчислень, для яких компʼютери власне і створювались. (Це здивувало всіх при першій взагалі комерціалізації компʼютерів, коли вони вийшли на початку 1950-х з воєнної секретности: замовлення під фінансові задачі були разів в 5 більше, ніж під математичну фізику.) Ми всі чули, що є такий COBOL і на ньому досі стоїть світова фінансова система, але це у першу чергу стиль його роботи з числами, який не дуже чисто (і ефективно) відтворюється сучасними засобами, що ростуть з інших доменів.

Як там робити, де і наскільки рухома крапка придатна — розберемось.

Цим список підходів не обмежується. В деяких випадках краще також використовувати раціональні числа, тобто дроби. Це вже за межами цієї всеї розповіді.

Загальне представлення і перші наслідки

Всі такі числа представлені єдиноманітно у всіх реалізаціях, незважаючи на безліченну кількість інших деталей:

± мантиса × (основа ↑ порядок)

крім спеціальних варіантів (нескінченність, не-число), про які треба вести мову окремо.

При компʼютерних обмеженнях використовують позначення з «e» («E») між мантисою і порядком, наприклад, 6.02e23 замість 6.02×10↑23. (Це не єдиний варіант, далі буде, проте головний.) Використовуватимемо таке позначення далі.

Але вже тут починаються проблеми деталей означень термінолоґії, що викликані довгим і частково незалежним розвитком.

Простіше всього зі знаком: він всюди знак (sign). Формула в показаному вигляді дає додатні (точніше, невідʼємні, якщо врахувати нуль; в більшості представлень є нуль з обома знаками), і відʼємні (недодатні) числа. Можна було б писати його більш математично, просто як «±» в загальному випадку, «+» або «-» у конкретного числа... але ми рухаємось до компʼютерного представлення, тому одразу використовуватимемо знак як однобітове значення: 0 — плюс, 1 — мінус; формулою це пишуть як (-1)↑s. В компʼютерному представленні знак — майже завжди 1 біт: якщо він 0, то число невідʼємне (0 або додатнє), якщо 1 — то недодатнє (0 або відʼємне).

Основа (анґл. «base», рос. «основание»; проте в ISO стандартах LIA, тобто ISO/IEC 10967, «radix») майже безальтернативно дорівнює 10 у не-IT використанні, і 2 у більшости застосувань у сучасних компʼютерах, але ви це можете і не побачити, якщо не зазирнете у деталі. Іноді бувають більш специфічні варіанти, як 8, 16; по більшости це спадщина старих реалізацій, частина характеристик така ж, як у основи 2, але через гіршу стабільність похибки зараз не рекомендується.

(Чому «майже» безальтернативно 10 за межами IT? А ви подивіться, як ми рахуємо час. В добі 24 години, в годині 60 хвилин, в хвилині 60 секунд. Раніше і секунди так ділили на терції, терції — на кварти, 60 кожної в більшій одиниці, і так далі... але зараз не ділять. Це мішана система з основою, що не дорівнює 10. Місяці з нерівним числом днів вже і не згадую. Ми звикли до цього, як Журден, і якщо не звернути увагу, не помічаємо специфіку.)

Мантиса — у деяких старих джерелах тут використовували термін «коефіцієнт». Анґлійською теж можна знайти, де він coefficient. Несподівано, в General Decimal Arithmetic Specification (про неї далі буде) теж coefficient. В чому проблема з терміном mantissa? Іноді його використовують, фактично, як лоґаріфм від того, що є мантисою в формулі: «дрібна частина лоґарифму числа». Наприклад, у числа Авоґадро за цим визначенням мантиса дорівнює 0.77975 (десятковий лоґарифм від 6.022...)

Мабуть, це стало причиною, чому при розробці IEEE754 (про нього буде далі дуже подробно, поки що просто згадуємо) створили окремий штучний термін significand. Проте він не став єдиним. Наприклад, в стандартах ASN.1 від ITU, в CBOR від IETF, цю частину звуть «mantissa» — і це при тому, що CBOR навіть замислили розробляти, коли IEEE754 був вже безальтернативним. Орієнтуватимусь на варіант «mantissa» навіть при «significand» в IEEE754, використовуючи останній переважно для того, як ця частина записана в бітовому представленні.

(Точна цитата з Changelog від General Decimal Arithmetic Specification: «The significand of a number has been renamed from „integer“ to „coefficient“, to remove possible ambiguities.» Нема проблем зрозуміти, що саме тут сказано, і до якогось ступеню навіть чому саме такими словами, але від формулювання дах готовий улетіти кудись в компанії качок у теплі краї...)

Часто мантиса нормалізована, тобто заведена у якісь межі. Є три основні стилі нормалізації:

• Найбільш популярний зараз: 1 ⩽ M < b, де b — основа. Приклад: 6.02 × 10↑23.
• Дуже популярний і місцями основний в раніші роки: 1/b ⩽ M < 1. Приклад: 0.602 × 10↑24.
• Специфічний для десяткової арифметики і для деяких зовнішніх представлень: мантиса — ціле число, що не ділиться на основу. Приклад: 602 × 10↑21 (якщо у нас нема більше значущих цифр).

У визначенні порядку їм відповідають left units view, fraction view и right units view; не перекладатиму ці терміни, хоча, може, хтось запропонує переклад. Про них буде далі.

Для нормалізації мантиси потрібна відповідна корекція порядку: 602 × 10↑21 = 6.02 × 10↑23 = 0.602 × 10↑24.

Порядок — анґлійською переважно зараз «exponent». Але можна зустріти і варіант «characteristic». Ще є варіант «scale», як у BigDecimal класі Java (причому з інверсним знаком), і в типах NUMERIC і DECIMAL в SQL. Слово «characteristic» дуже довге; мабуть, тому «exponent» поступово перемагає. Але його не треба плутати з exponent у загальному випадку піднесення чогось в степінь. Ми використовуватимемо слова «порядок» і «exponent».

Тепер, коли читач остаточно заплутався, можна перейти к наступній частині. :irony:

Математичні дійсні числа, природні, ті, що в компʼютері, і ті, що у нашому мозку

Для чистої математики не може бути обмежень на те, чому рівне дійсне число... поки воно дійсне (не ціле, не раціональне; і не комплексне, бо це інша казка). Множина дійсних чисел підкорюється ряду аксіом і наслідків з них. Наприклад, між двома різними дійсними числами завжди можна вибрати ще одно, яке не дорівнює обом — а, як наслідок, нескінченну кількість таких чисел. А якщо ще подумати — не просто нескінченну, а контінуум. (Хто не знає теорію множин — пропустіть останнє, для нашої розмови непринципово.) І порядок числа необмежений в обидва боки, навіть якщо для таких чисел немає фізичного еквіваленту.

Для фізики вже вступають в дію обмеження квантової механіки, але не всюди (у макро-світі вона дає менший вплив, ніж iнші похибки — власне, саме так і визначають межі, де треба починати думати про квантову механіку).

Але в компʼютер ви так просто не помістите числа довільної точности: памʼяти не вистачить. І цього ніхто не робить ще й з точки зору виробности: треба, щоб результат був отриманий за скінченний час:) Тому обидва показники: точність, у кількости цифр вибраної основи, яка може максимально вміщуватись в значеннях, і діапазон порядку — обмежені у конкретній реалізації.

А це означає, що замість точних значень ми завжди працюватимемо з їх наближеннями, що вибрані з якоїсь множини доступних значень. Замість неперервної шкали маємо набір точок, і для кожного значення маємо обирати одну з близьких до нього точок. Шкала неоднорідна, ближче до нуля точки стоять купніше, ніж даліше від нуля. Наприклад, десяткова при 3 цифрах точности матиме: 9.95, 9.96, 9.97, 9.98, 9.99, 10.0 (був крок 0.01), а далі 10.1, 10.2, 10.3, 10.4... (крок вже 0.1). Ну і сам набір цих точок обмежений (по модулю) як знизу (навколо нуля), так і зверху (у напрямку нескінченности).

У частині 2 будуть «неприємні» наслідки з цього.

Використання в мовах проґрамування, основи

Ми маємо фактично дві ґрупи мов: одні наслідують максимум синтаксису від C (включаючи те, що давно визнано помилками), інші ідуть своїм шляхом.

Якщо це C-like (C, C++, Java, C#...), то у них будуть слова float для одинарної точности і double для подвійної. Примітимо, не single і double, хоча сам тип в дотнеті зветься Single. На зараз, якщо ви не використовуєте щось дуже специфічне, вони будуть 32- і 64-бітними двійковими типами IEEE754. (Ще може бути long double, з ним складніше, залежить від платформи.)

Далі все прямолінійно — імпортуєте дані (або берете константи), додавання, множення, віднімання, ділення, ітд., як звично, тільки символи можуть бути інші, і за ґрупуванням треба слідкувати, і вставляти дужки де треба. Одного прикладу, сподіваюсь, досить:

// returns root amount; 0, 1, or 2 are exact amount;
// 3 means infinity of them, no x1, x2 assigned.
// Does not handle overflows, underflows, precision loss.
int
the_most_naive_square_equation_solver(
    double a, double b, double c,
    double *x1, double *x2)
{
  if (a == 0) {
    if (b == 0) { return c == 0 ? 3 : 0; }
    *x1 = *x2 = -c/b; return 1;
  }
  double d = b*b - 4*a*c;
  if (d < 0) { return 0; }
  if (d == 0) { *x1 = *x2 = -b/(2*a); return 1; }
  d = sqrt(d);
  *x1 = (-b+d)/(2*a); *x2 = (-b-d)/(2*a); return 2;
}

В Rust, Zig імена двох основних типів — f32, f64. В Go — float32, float64. Новітні мови схиляються до такого явного іменування. Операції з даними, в цілому, ті ж самі.

В Fortran одинарна точність це REAL, подвійна — DOUBLE PRECISION; але є можливість писати в стилі REAL(8), REAL*8, залежно від стандарту і діалекту. Є і більш дивні форми, наприклад, «REAL, KIND(1.0D0) :: X» визначає X як змінну подвійної точности.

В JavaScript є один number, він 64-бітний, а фактично це єдиний вбудований чисельний тип (деякі операції з ним, які сприймають як 32-бітні цілі, за межами цеї теми). В Python (без numpy) є один float, він 64-бітний. (Numpy вже знає float16, float32, float64; місцями і float128, але з ним проблеми, чи дійсно він 128.) Оскільки веб у всіх його видах це більше половини всього проґрамування, одинарна точність в цій субгалузі майже не використовується. Єдиність числового типу, і те, як з ним працюють, реґулярно створюють ситуації, коли на виклик таксі вам приїжджає [Object object] номер undefined через NaN хвилин.

В стандарті SQL є REAL, який «implementation-defined» одинарна точність, DOUBLE PRECISION, який відповідно «подвійна» точність (не ґарантовано), і FLOAT, якому можна задавати точність через параметр, у бітах; але є і окремий набір NUMERIC і DECIMAL, яким точність задається у десяткових цифрах (їх краще вважати типами з фіксованою крапкою).

Ви познайомились з нашим звичним безладом (точніше, з одним його боком), продовжуємо нашу екскурсію.

Типи потрібні не тільки для представлення даних, а ще й для операцій з ними. В стандартному комплекті (орієнтуюсь для зразка на glibc) є:

1) Арифметика — додавання, віднімання, множення, ділення.

У переважній більшости мов, які зараз в ходу (COBOL не враховуємо, з його «ADD X TO Y GIVING Z», і то, з якогось дуже давного часу є COMPUTE з більш звичним синтаксисом), для складання треба написати a+b, для ділення — a/b... Стандартні операції і, по більшости, очікувані їх результати.

Специфіка де-де: ціле ділення (з залишком) і плаваюче можуть позначатись різними знаками (Python, Pascal/Modula/Ada) або розділятись згідно типам арґументів (C, C++, Java, C#...)

2) Степені, корені, лоґарифми. Вже не знаками операцій, а функціями з іменами, як і більшість того, що далі.

3) Триґонометрія. sin, cos...

4) Більш рідкі математичні функції — навіть такі, як ґамма-функція, функція помилки (error function — erf), інші дива. Нема в IEEE754, але є в стандарті C.

5) Порівняння.

Операції класів, що були до цього місця, пишуться явно. А ось далі...

6) Конверсія між внутрішнім і зовнішним представленням («ввід-вивід», «імпорт-експорт»), у різних форматах з різною точністю. Уникатимемо термінів "ввід"/"вивід«, бо вони приймають на увазі взаємодію з файлом чи зовнішнім пристроєм; імпорт і експорт — більш нейтральні.

Імпорт дуже часто схований за різними там parseNumber() і scanf(), або навіть +x (JavaScript). Експорт — printf, toString(), ""+x(JavaScript), форматні рядки (f"{var1} {var2}" в Python, аналоги в інших). Або навіть взагалі не видно, що експорт чи імпорт має місце — обробка проходить автоматично. Наприклад, в шеллах (sh, bash), Tcl, єдиний тип даних — рядок, розрізняється тільки його обробка. І тут можемо за рахунок неочікуваної конверсії (або неочікуваного режиму конверсії) щось утратити.

7) Конверсія типів (між цілим і плаваючим, між різними розмірами). Сюди ж — цільові округлення.

Це теж конверсія, і вона може проходити непомітно. Написали f = i; — ось вам і конверсія.

8) Допоміжні операції — класифікація, масштабування, розділення на мантису і порядок або зборку в число, керування контекстом.

Взагалі, поки тонких ефектів нема (у 99.9+% випадків того, з чим зустрічається пересічний проґраміст), все виглядає просто і стабільно. І знову — це треба викарбувати на камені.

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

Математичні дійсні числа, природні, ті, що в компʼютері, і ті, що у нашому мозку

В частині 1 говорилось про обмеження точности чисел у реальних представленнях. Які наслідки з цього?

Звичайний формат одинарної точности IEEE754 дозволяє представити числа від 8388608 (2↑23) до 16777216 (2↑24) з кроком 1. Тобто за 8388608 іде 8388609. Ви не зможете в ньому представити 8388608.5 чи 9876543.746. Після 16777216 крок вже 2. Тобто 16777216 і 16777218 представити можна, а 16777217 — ні.

Да, це одинарна точність. А що в подвійній? Перше натуральне число, яке не можна представити, дорівнює 2↑53+1, тобто 9007199254740993. 16 цифр. (Хтось тут вже побачив, що номера платіжних карток можна представити як число в JavaScript чи PHP. Якщо починаються з 4 або 5, як майже всі у нас, і вміщаються в 16 цифр... як казав один мій колеґа, «я ще розумію, як можна так писати, але як можна таке рекомендувати???» Ні, я не рекомендую. Бо не помітите, як щось зіпсується, і воно перестане працювати.)

Точність... ми всі (ну добре, майже всі) думаємо в десятковій системі. Наприклад, «тут вірні 5 цифр після крапки». Традиційні посібники з метролоґії, фізичних експериментів розраховані на це, навіть на шкільних лабораторних роботах навчають такому стилю: плюс-мінус не, наприклад, 0.003125, а 0.001, якщо можете, або вже 0.01. Ось берете і пишете: маса вантажу 4.91±0.01 кґ. Багато законодавчих актів і галузевих реґуляцій сформульовано по тим же принципам. (Наприклад, маса яйця 1-ї катеґорії від 53 до 62.9 ґрам, а 0-ї від 63 до 72.9 — очевидно, що її треба рахувати з точністю до 0.1 ґрама, а округлювати до цієї точности за «фінансовим» стандартом.) Це настільки вʼїлось, що стандарти C і C++ дають константи для цеї точности для наявних типів з рухомою крапкою. (Здавалося б, мови близькі до нижнього рівня... а туди ж.) Що ми можемо сказати про точність у конкретних типів? Знову беремо IEEE754 (хоча специфіка загальна для всіх представлень з основами 2, 8, 16).

Тут починаються особливості того, що представлення двійкове, а не десяткове.

Число з одинарною точністю (binary32) IEEE754 (як воно представлено — буде далі) може передати (якщо результат скінченний і нормалізований) довільне десяткове число до 6 значущих цифр при довільному порядку. Є випадки, коли число з 7 цифр при імпорті (конверсії в двійкове) і експорті (назад в десяткове) не може бути збережено (два таких різних числа дають однаковий результат): наприклад, 8.589973e9 і 8.589974e9 дають одне і те ж (при зворотньому експорті буде друге з них).

Але: щоб передати довільне число IEEE754 з одинарною точністю в десятковому вигляді, потребується до 9 (включно) десяткових цифр. Є випадки, коли при меншій кількости цифр початкове число не буде збережено: 100000021.0; 0.102000006; 10.0060005... (проґрамний пошук починався з чисел, як 0.1, тому такі результати. можете почати з будь-якого іншого).

6 і 9. Різниця кричуща. Хоча і сам факт, що ці два значення не рівні, вже може збити з пантелику непідготовленого.

Для подвійної точности тут буде 15 і 17. (Дивіться вище про номери карток.) Різниця дещо менша. Але 2 це теж не те, чого бажалось в ідеалі. Насправді, менше 2 її зробити неможливо.

Краще не буде. Далі буде розписано. Поки що просто відмічаємо собі: тут мешкають дракони, і вони дійно можуть вкусити.

Але... люди продовжують думати в десятковій системі. І це принципово.

Математика, що в ній непросто.

Колись знаменитий FMM не втратив за майже 50 років своєї актуальности тим, що дає сконцентрований список типових проблем, які виникають в математичних розрахунках при переході від ідеального світу (без похибок, округлень і лімітів представлень) до реальних обчислень:

1. «Некоректність» самої задачі: мале відхилення вхідних даних викликає довільно великі (математично, такі, що не можна обмежити якимось корисним множником) зміни вихідних даних.

В FMM наведені приклади таких задач, починаючи з поліноміального рівняння одної змінної (x-1)*(x-2)*...*(x-20)=0, де всі корені дійсні і рівні 1, 2, ..., 20; зсунули коефіцієнт при x↑19 на 2↑-23, тепер десять коренів вже комплексні.

Ми також можемо зіслатись на проґнозування погоди: ця задача некоректна, мінімальне відхилення в одній точці може розповсюдитись на всю планету, точне передбачення навіть до факту «буде дощ чи ні» більш ніж на ≈3 дні неможливе. Відкриття цього факту в 1960-х було шоком.

Штучний, але красивий приклад — Mullerʼs Recurrence. Якщо проміжний результат перевищить правильне значення 5 хоч на хвильку, обчислення стрибне до наступного стабільного значення 100. А таке перевищення виникає у більшости випадків через похибку. Примітьте спеціальний характер прикладу: проблема виникає при будь-якій точності! (може, повільніше, але все одно вибухає) — бо і нові відхилення виконуються в тій же точності.

Деякі задачі можуть бути «стабілізовані» зміною постановки, так звані «методи розвʼязку некоректних задач» розробляються вже пару століть. Сподіваюсь, що ваші задачі не мають таких проблем (ви не вираховуєте погоду;))

2. Переповнення або уход в нуль проміжних результатів.

Простий приклад... проґрама «банального» (всі вчимо в школі) квадратного рівняння. x↑2-2*x+1=0 розвʼязати тривіально? А якщо 1e300*x↑2-2e300*x+1e300=0? Привіт, при обчисленні в стандартному double, в b↑2-4*a*c спочатку нескінченність зменшуваного і відʼємника, а потім не-число (бо нескінченність мінус нескінченність). Тепер спочатку нормалізуємо коефіцієнти, у сенсі, до мінімального за модулем спільного порядку, потім вже починаємо щось кудись возводити. Кажете, штучний приклад? Якщо на вході космічні відстані, які з меґапарсеків перекладені у метри, ані разу не штучний.

3. Повна або суттєва утрата точности при складанні великих (значно більших за результат) взаємно компенсованих доданків.

Є приклад навіть для «того ж» квадратного рівняння, x↑2-100000000*x+1=0. Якщо вирахувати менший корінь напряму через повну формулу, а точність одинарна, отримуємо 0. Ну, тут я не знайшов «реального» виправдання такому прикладу. Але і не треба повної втрати точности, щоб побачити проблему, достатньо 1-2 найменших розрядів. Тому вираховуємо спочатку більший за модулем корінь, а менший — не довгою формулою, а діленнєм (x2 = c/(a*x1)).

Або обчислення степеню при значному відʼємному показнику. Краще вирахувати 1/exp(20), ніж exp(-20), через звичайний ряд-сумму, бо поки великі складники почнуть впадати, похибка в менших розрядах убʼє точність результату. Ну, це перший рівень складности; щоб отримати точний результат згідно вимог IEEE754, таких простих прийомів недостатньо.

З новішіх аналоґів FMM є, наприклад, Accuracy and Stability of Numerical Algorithms (2002), але його складніше знайти; FMM (у російському перекладі) доступний «вільно».

Сподіваюсь, більшости читачів цеї статті не знадобиться глубоко влазити в тематику саме математичних проблем. Але згадати було потрібно — бо треба памʼятати, що таке в принципі можливо, і своєчасно відреагувати.

Граблі загальні, але дитячі...

В десятковій системі числення можна точно представити будь-який дріб зі знаменником, що ділиться на степені 2 і 5 (що значить і 10), тому, наприклад, 0.1 там представлено точно. В двійковій системі це урізається тільки до степеню числа 2. 0.1 в двійковій системі виглядає як 0.0001100110011..., або 0.0(0011) в тих же позначеннях, в яких в десятковій 1∕3 пишеться як 0.(3). Ніякий скінченний позиційний дріб це не представить, і питання в тому, яку довжину ми обріжемо і як округлимо.

Робимо простий текст (на Python, просто для швидкости):

cnt = 0
i = 0.0
while i <= 20.0:
    cnt += 1
    i += 0.1
print(cnt)

Можна повторити на будь-якій мові, вибравши подвійную точність, тільки не вмикайте fast-math і тому подібні опції небезпечних оптимізацій.

Який результат можна наївно очікувати? Не забуваємо про ефект плюс-мінус одиниці (наприклад, ряд значень 0, 1, 2, ..., 9, 10 — це 11 значень, а не 10). 201? Ні, воно надрукує 200. Де ще одна ітерація, хто її зʼїв? Дивуємось, додаємо print(i) після циклу. Воно показує: 20.000000000000014. Звісно, це більше 20.0, тому останнє порівняння «<= 20.0» не виконалось і одну ітерацію було пропущено.

Той же ефект для кроку 0.01: 2000 зроблених ітерацій замість 2001.

А якщо ліміт з 20.0 замінити на 10.0, ефект кудись зникає, маємо 101 ітерацію при кроку 0.1 (і 1001 для 0.01). А фінальне print(i) друкує... 10.09999999999998 (10.1 без малого). Додамо друк попереднього значення суми. Для кроку 0.1 і ліміту 10.0 там буде 9.99999999999998.

Тобто додавання в діапазоні до 10 дає легку недостачу значення (і тоді умова виходу коректно спрацьовує), а до 20 — вже перебільшення (і тоді — ні).

Чому такі ефекти? Бо 0.1 не можна точно представити у двійковому вигляді. Якби у нас була необмежена точність, ми б його представили як 0.0(0011)₂, де дужки показують повтор — таким же чином, яким 1/3 в десятковій позначається 0.(3). Але у нас точність обмежена. 0.1 в подвійній точності передається як 1.1001100110011001100110011001100110011001100110011010₂ домножене на 2↑-4 (тобто 1/16). Примітимо, що хвіст тут округлений (інакше б дві останні цифри були 01).

Всі числа неточні — тому всі суми неточні. А далі — як карта ляже. Може округлити сумму в плюс, може в мінус. Сусідні округлення можуть скомпенсувати одне одного, а можуть, навпаки, скластись зі збільшенням ефекту. Хто хоче, може сам погратись с різними кроками і лімітами, спробувати також з 32-бітним типом замість 64-бітного, як було у Python... У кожному конкретному випадку можна зрозуміти, що діється. Але передбачити так, щоб обійти ефекти — надто дорого і тому непродуктивно. Як універсальний принцип, краще за все запамʼятати: таку ітерацію по float робити не можна, буде накопичення мікропохибок.

(Або, у особливих випадках і ретельно контролюючи наслідки, замінити порівняння, наприклад, на <= 20.05 замість 20.0. Треба при цьому довести, що не буде накопичення похибки. Все одно краще так не робити.)

Після цього вже не буде настільки дивувати цей відомий мем:

Спитаймо Python про внутрішній формат такого числа порівняно з тим, що він сам вважає 0.3:

>>> binascii.hexlify(struct.pack('>d', 0.3))
b'3fd3333333333333'
>>> 0.1 * 3
0.30000000000000004
>>> binascii.hexlify(struct.pack('>d', 0.1 * 3))
b'3fd3333333333334'

Як розуміти ці результати більш усвідомлено — буде в частині 3. Поки що достатньо розуміння, що це представлення деякого числа в шістнадцятковому вигляді, і бачити, що вони відрізняються одною найменшою цифрою. Різниця тут в один Unit of least precision (ULP — цей акронім часто зустрічається в тематиці), тобто найменшу цифру з представлених.

І це все тільки верхушка безодні. Трошки-трошечки копнути і спробувати повернути початкове...

>>> binascii.hexlify(struct.pack('>d', 0.3/3))
b'3fb9999999999999'
>>> 0.3/3
0.09999999999999999
>>> binascii.hexlify(struct.pack('>d', 0.1))
b'3fb999999999999a'
>>> binascii.hexlify(struct.pack('>d', (0.1*3)/3))
b'3fb999999999999b'
>>> (0.1*3)/3
0.10000000000000002

Три числа, між кожною парою сусідніх 1 ULP, але вони самі неочікувані. Ми домножили 0.1 на 3, розділили назад, отримали не те, що 0.1, з якого починали. А 0.1, в свою чергу, не дорівнює ні тому, що 0.3/3 (яке теж інше, ніж (0.1*3)/3). «Ось і живіть з цим» (™).

(Насправді, тут і показ дещо спрощений. Щоб результат 0.1×3 був показаний як 0.30000000000000004, а само 0.3 не було показано як 0.29999999999999999, треба спеціально скоротити представлення 0.3; і 0.1 з такою точністю це не 0.1, а 0.10000000000000001. Якщо ви бачите компактні форми 0.1, 0.3, вони вже скорочені спеціальним алгоритмом без втрати точности. JavaScript, Python це роблять зараз, але багато хто не робить. Це обговоримо в частині про зовнішні представлення. Поки що приймемо це все як належне і радіємо, що нам показують цю різницю.)

І в іншу сторону:

>>> 0.3*3
0.8999999999999999
>>> 0.1*9
0.9
>>> (0.1*3)*3
0.9000000000000001

Останнє ще й показує, що операція множення тут не ассоціативна.

Проте іноді очікувана відповідність відроджується сама:

>>> 0.3*3 - 0.2
0.7

А також згадаймо цю історію:

> The software error of a MIM-104 Patriot caused its system clock to drift by one third of a second over a period of one hundred hours — resulting in failure to locate and intercept an incoming Iraqi Al Hussein missile, which then struck Dharan barracks, Saudi Arabia (February 25, 1991), killing 28 Americans.

(цитується по вікі).

Накопичуєте без контроля і корекції — отримуєте результат.

Колись так мій колеґа нарвався з циклом по float саме так, як я показав в першому прикладі з <=20.0, і дуже дивувався... а я теж був шокований тим, що він цього не знав, і зміг тільки сказати «це ж в книгах на перших сторінках!» Колеґа образився... але, сподіваюсь, більше не на мене (аберація, згоден — я вчився на прикладну математику, а він на дещо ближче до електроніки), а на саму ситуацію, коли те, що здавалось простим і очевидним, раптом вибиває ґрунт з-під ніг. Дійсно, рухома крапка — вона така...

Використання в мовах проґрамування, проблеми.

В першій частині були принципи використання, у позитивному ключі — тобто воно є і з ним можна працювати. Переходимо до зворотнього боку.

1) Арифметика — додавання, віднімання, множення, ділення.

Результати (a+b)+c і a+(b+c), або a*(b+c) і a*b+a*c, можуть бути нерівні і давати зміну результату — звісно, тут знаками + - * / і так далі позначені операції, як вони виконуються в компʼютері, а не їх теоретично ідеальні аналоґи. Це зветься, що операції, як додавання, віднімання, множення, ділення неассоціативні (віднімання, ділення — з відповідною корекцією на зворотні операції), множення недистрибутивно. (В математичних описах, переважно, розділяють при цьому знаки: «+» — ідеальне, а «⊕» — реальне у проґрамі, і так далі для всіх операцій; тут відкладатиму їх до останнього.) Наприклад, a*b і a*c обидва настільки малі, що вироджуються в 0, аa*(b+c) має ненульовий результат. Або вище був приклад, як a*(b*c) ≠ (a*b)*c, коли a = 0.1, b = c = 3.0. Тому компілятори без прямої вказівки проґраміста не роблять такі заміни, на відміну від цілих чисел. З GCC треба додати `-ffast-math`, щоб дозволити такі ігри з виразами. В Java до якоїсь версії треба було додавати `strictfp`, щоб такого не траплялось. В сучасних версіях цього вже не треба, вирішили, що краще мати strictfp завжди.

2) Порівняння. Це вже специфіка не завжди, а IEEE754, але у нас зараз альтернативи немає: кожне з порівнянь {1 < NaN}, {1 == NaN}, {1 > NaN} є хибним. Про це буде окрема частина.

Операції класів, що були до цього місця, пишуться явно. А ось далі...

3) Імпорт-експорт, тобто конверсія між внутрішнім і зовнішним представленням, по більшости десятковим.

Рядки «3 » і «3» (в першому в кінці простий пробіл) це одне і те ж — представлення числа 3, чи ні? Bash, Python, JS скажуть при спробі імпорту, що да. А якщо після цифри іде U+2002 (EN SPACE)? Bash вже не згодний, що це число; Python і JS ще згодні. А якщо після цифри іде U+200B (ZERO WIDTH SPACE)? Python вже не погодиться, а для JS це все ще число. А «010» це 10 чи 8? Хто ґарантує читання як float, щоб не спрацювало правило сішного леґасі першого післяпотопʼя «префікс 0 означає вісімкове число»? Чи треба розуміти «123z» як 123, або ж вважати імпорт невдалим? Чи буде читатись «123_456» так же, як «123456»? Локалі: «123.456» це 123 цілих 456 тисячних (C, en_US, багато інших) чи 123 тисячі 456 цілих (da_DK)?

4) Конверсія типів (між цілим і плаваючим, між різними розмірами). Сюди ж — цільові округлення.

Це теж конверсія, і вона може проходити непомітно, а іноді давати неприємні ефекти. В проґрамах на C/C++ на більшости актуальних платформ (і не embedded) int — 32 біти. Точність 32-бітного IEEE float — 24 біти. Конверсія автоматична, якщо в операції другий арґумент float, утрата точности проходить непомітно. Ось Брюс Доусон лаявся на це. (Взагалі, рекомендую його блоґ тим, хто любить «системне» проґрамування. На жаль, нових технічних постінґів не буде, але і попередніх вистачить на багато днів.) Але те ж саме маємо, якщо з 64-бітного long зробити double: 53 біти замість 64 — якось замало. Присвоїти значення типу double в змінну типу float — теж мовчки зʼїсть. На жаль, в цих мовах скасовувати неявну конверсію чи обкладати її складними допоміжними умовами вже запізно; але в нових мовах (у першу чергу згадався Rust) це роблять.

(І ще більше тонкощів: ось в формулі коренів квадратного рівняння є 4*a*c і 2*a. Якщо ми по-тупому вимагатимемо явної конверсії всіх цілих в плаваючі, отримаємо неможливість так писати: треба або 2.0*a, або double(2)*a... Тому в мову краще додати правило: константа, як 2 або 4, до якогось моменту нетипізована, і може представляти і ціле, і плаваюче. Такі правила є зараз в багатьох. І це показує складність, яка стоїть за зовнішньо простим кодом.)

В наступній частині: головний поточний стандарт IEEE754, що, нащо, чим гарний і так далі.

Підписуйтеся на Telegram-канал «DOU #tech», щоб не пропустити нові технічні статті