Плаваюча (рухома) крапка, частина 4: Сучасний стандарт — тонкі деталі і проблеми

Читайте також:
— Плаваюча (рухома) крапка. Частини 1-2: ввідна, перші граблі. Там же каталог частин циклу.
— Плаваюча (рухома) крапка. Частина 3: Сучасний стандарт і все навколо — основи
— Плаваюча (рухома) крапка. Частина 5: Текстовий експорт-імпорт
— Плаваюча (рухома) крапка. Частина 6: Порівняння
— Плаваюча (рухома) крапка. Частина 7: десяткова рухома крапка

Це дуже мала частина, але краще виглядало відокремити її.

В частині 3 була база даної підтеми. Перейдемо до тонких і спірних моментів (не повʼязаних тут з взаємодією з людиною).

Сінус, такий собі сінус... і не тільки він. Безцільна точність?

Почнемо з такого дивного прикладу: функція sin() має період рівний 2×π ≈ 6.28... Сітка значень, які можуть бути представлені як число з одинарної точністю, така, що в діапазоні від 2↑26 = 67108864 до 2↑27 = 134217728 представляється із кроком 8 (67108864, 67108872, 67108880...), а те, що між ними, не можна представити. Тобто, якщо таке число, наприклад, 67108872 подане на вхід sin(), то ми його знаємо з плюс-мінус 4 (повний розмах похибки дорівнює 8), його точність вже втрачена за рахунок попереднього округлення, і ми могли отримати будь-який результат в діапазоні значень sin(), від —1 до 1. Чи є взагалі сенс вираховувати значення sin() від такого арґумента?

Я б сказав — ні, сенсу вже нема, втрачено все, що могли; якщо sin() віддаватиме для такого арґументу завжди 0, це нічого реального не спотворить. Але стандарт вимагає точного вираховування для того арґументу, який був переданий. Це означає, що арґумент буде спочатку нормалізований, з повною необхідною точністю, до вхідного діапазону [0; π∕2], і потім вже «цивілізовано» вирахований результат з такого арґументу. Що це означає для реалізатора?

При значенні арґументу аж до 2↑128-2↑104 = 3.40282347e38 (максимальне для одинарної точности), тобто повних 128 біт цілого значення, і необхідних 24 дробової після нормалізації, максимальна довжина обчислень, що може знадобитись для такої нормалізації, дорівнює 152 бітам. Для подвійної точности це вже буде 1077 біт. (Для четверної точности, тобто в 128-бітному форматі — 16 497!) До цих цифр треба буде ще додати 2-3 біти на коректне округлення у процесі (мінімум згідно правилам, описаним в попередній частині) і ще, можливо, на проблему наступного підрозділу, але вже названі цифри лякають.

Як казали в одному анекдоті, «Так, жах. Але не „Жах, жах!“» У більшости випадків цього не треба. Фактичний код (наприклад, тут) розділяє прості і складні випадки, і підіймає таку машинерію тільки для великого (за модулем) арґумента. Але особисто мені все одно не до кінця зрозуміло, для чого взагалі такі випадки обробляти з максимумом точности. Ну, вже стандартизували, алґоритми є, не видаляти ж...

Table-maker dilemma і точність трансцедентних функцій

Ще одна проблема вимоги «нескіченно точний результат, округлений згідно режиму», щільно повʼязана з попередньою, виглядає так. Абстраґуємось від стандарту, і ось нехай нам треба округлити, для прикладу, до цілих, результат деякої хитрої функції. Округлення half-to-even, і 2.5 округлиться до 2, а якесь 2.500001 — до 3. Ми рахуємо значення, і воно поки що рівно 2.5, 2.50, 2.500, 2.5000..., ґенеруємо ланцюжок нулів, і ми не знаємо, буде воно рівно 2.5, чи можна зупинитись, або треба продовжувати, поки не побачимо, що там таке є ненульовий хвіст (і тоді результат має бути 3). Коли зупинятись? Для будь-якої штучно вибраної точности знайдуться функції і арґументи, для яких цеї точности не вистачить.

Приклад, звісно, грубий і штучний до неподобства, реальні будуть значно тонкіші. AI мені знайшов, наприклад, cos(5992555.0) в double. Або з одної статті:

• log10(2.60575359533670695e129) = 129.41593334566474 (найточніше) чи 129.4159333456647 (деякі версії Microsoft C рантайма)?
• expm1(-1.31267823646623444e-07) = —1.3126781503100307e-07 (найточніше) чи —1.3126781503100305e-07 (бачу на своєму Linux)?

Або ось порівняння різних реалізацій на прикладі sin(). (Тут одразу і лікування.)

Це дає проблему коректности результату з точністю, у найкращому випадку, до одного знаку найменшої представленої точности, unit of least precision (ULP): плюс-мінус один такий ULP — це те, наскільки можуть відрізнятись різні реалізації одної функції для одного і того ж арґументу... а може, і більше. (Але це саме в термінах плюс-мінус одної цифри низького порядку. Не влізатимемо подробно, але, якщо це різниця між, десятковим прикладом, 1.99999 і 2.00000, формально ми можемо не вгадати жодної точної цифри! Це вже проблема формального визначення, що ми вважаємо точністю.)

Intel пише в своєму SDM: «The P6 family and Pentium processors’ results will have a worst case error of less than 1 ulp when rounding to the nearest-even and less than 1.5 ulps when rounding in other modes.» Фактично це має бути і для сучасного FPU (зватимемо FPU в x86 процесорах як x87 FPU, згідно його походженню). Для SSE ця норма напряму непридатна, він не підтримує степені, лоґарифми і триґонометрію, вони всі віддані на рівень софта. Але в софті чи в залізі, проблема завжди буде.

Те що в попередньому абзаці — для FPU, як вже зрозуміли, відхилення від стандарту. Але і стандарт тут пише тільки про «recommended» operations, хоча, «A conforming operation shall return results correctly rounded for the applicable rounding direction for all operands in its domain.» — тобто, якщо ви дозволяєте похибку в 1 ULP через згадану проблему, то це вже не conforming. Розгадка: стандарт тут спирається на факт наявности з 2001 року зразкової реалізації степеней, коренів, лоґарифмів і триґонометрії — IBM Accurate Portable Mathlib. Після цього зникли найменші сумніви, що це може бути досягнено... але, велике «але»:
1) тільки для обмеженого кола функцій,
2) реалізації вивірені тільки для фіксованого набору точностей IEEE754,
3) іноді надто дорого (і знову згадуємо embedded).
Звісно, легко додати вимогу, коли ти обмежуєш, до чого вона застосовується, тим, що вже розвʼязано... (трошки сарказму)

На практиці це може бути помʼякшено (якщо критично мати швидке виконання операції), але при дотриманні іншої вимоги: похибки дозволені за «best effort», поки функція залишається монотонною на тих же інтервалах арґументів, як і її «ідеальний» аналог (бо порушення монотонности створює значно більші проблеми!) Наприклад, лоґарифм це монотонно зростаюча функція; тоді якщо A ⩾ B, то необхідною вимогою є log(A) ⩾ log(B) для тої log(), яка реально реалізована. Це задовільнити значно простіше і дешево. X87 FPU ґарантує монотонність навіть при таких похибках, як у нього.

Всі такі реалізації занадто дорогі для деяких випадків: код впадає в проміжну точність аж до 768 біт (для double). З GNU libc, де спочатку вони були всі виконані кодом з IBM Accurate Portable Mathlib, ці випадки зараз видалені, починаючи з 2018 року... тобто вона, формально, відповідала, а потім перестала відповідати стандарту!😲 В вікі є подробиці зі списком реалізацій.

Дещо дивна назва «table-maker dilemma» натякає на те, що в практиці того часу, коли її сформулювали, проблема виражалась в тому, до якого рівня подробиць (скільки бітів в індексі) створювати таблиці. Думаю, зараз її назвали б якось інакше.

FPU не рекомендується для сучасного використання, в 64-бітних режимах треба особливі методи, щоб його застосувати; вся рухома крапка на x86 за замовчуванням виконується через SSE. Інші платформи загального призначення (в першу чергу думаємо про ARM AArch64, який зараз в 99+% мобільних пристроїв) теж покладають такі проблемні операції цілком на софт. Але в GPU нормально побачити апаратну реалізацію трансцедентних функцій... і там можна знову побачити схожі проблеми. Пошук підказує, що реалізація у NVidia може витрачати значно більше часу на особливі аргументи (це і до попереднього підпункту).

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

Висновок цього підпункту:
1) Якщо у вас не embedded, і немає потреби в ідеальній точності в межах доступної точности, вас це не хвилює.
2) Якщо embedded, то подумайте тричі, наскільки вам потрібні операції зі стандартної бібліотеки, часто простіші заміни достатні.
3) Якщо потрібна ідеальна точність... мої співчування, і... або міняйте бібліотеку на більш точнішу, або модель, щоб не була такою чутливою.

Не-число, чим саме воно не число

У випадках, типу квадратного корня з мінус одиниці, якщо не робиться миттєва ґенерація винятка з переходом на обробник, результатом є так зване «не-число» (not a number, NaN). Воно показує «ніякого значення тут (вже) нема». Те ж саме можна отримати, наприклад, операцією Infinity мінус Infinity.

Те ж саме вимушене спрощення стосується і значення «нескінченність», але воно показує, що ми приблизно знаємо, що трапився вихід за межі, але розрахунок був ще коректний.

Отже, у нас є нескінченність (Inf, Infinity) і не-число (NaN). Нескінченність одна на кожний знак (всього дві). Не-чисел більше; поступово склались два великих класи — сіґналізуючий і тихий NaN (signaling NaN ≡ sNaN; quiet NaN ≡ qNaN). В двійкових форматах IEEE754 найстарший біт significandʼу показує підтип: 0 — signaling, 1 — quiet. Що цікаво, в десяткових форматах строго навпаки (0 — quiet, 1 — signaling), що показує, що дехто виробив десятковий підстандарт до того, як фіналізували норми для двійкового (в 2008-му встановлено, а в 1985-му ще ні). Signaling як арґумент викликає ґенерацію помилки «invalid operation» при майже всіх операціях зі значеннями, quiet — не викликає.

У не-числа є так званий payload (як перекласти?), який не має значення для стандартних операцій. В ньому може бути що завгодно, хоча є рекомендації щодо canonical payload. У 1980-90-ті намагались виробити правила, як формувати цей payload, щоб передавати деталі невдалих операцій; але швидко здались, і десь після 2000 цю ідею вже не пропаґують. З іншого боку, ніщо не заважає користувачам помічати через нього різні початкові джерела «не-даних».

Але це була тільки перша частина... добре, ми розуміємо, що NaN + 1 = NaN, log2(NaN) = NaN, і все таке. А чи можна сказати, що NaN більше ніж 1? Чи менше?

І ось тут виступає на сцену цікаве рішення: NaN (довільного підтипу, з будь-яким payloadʼом) не порівнюється ні з чим — не менше, не більше і не рівне; повністю незалежне значення. (У тому числі з самим собою.) Це фундаментальне рішення впливає на багато аспектів; наслідки будуть розглянуті далі в розділі про порівняння. Тут же зупинимось.

Виконання простих операцій

Раніше було описано «ідеальне» виконання простих операцій за рахунок підвищеної точности проміжних результатів (до 3 біт або 2 десяткових цифр), RPSP (фон Неймана) до цеї підвищеної точности і фінального округлення. Здавалося б, все нормально? Але:

1. У деяких випадках кожне проміжне округлення тут щось помітно псує. Згадуємо знову цикл з кроком 0.1 до 20.

Нема зовсім загального вирішення цеї неприємности. Але є легкі допоміжні рішення, які можна застосувати у широкому класі випадків. На зараз є, мабуть, один такий: fused multiply-accumulate (FMA); вимагається в стандарті починаючи з версії 2008 року, є в сучасному X86, є в сучасних ARM, RISC-V, багато ще де... В ньому операція виду a+b*c виконується так, що вимога «так само як і з нескінченною точністю» розповсюджується і на це проміжне b*c, а округлення застосовується тільки після фінального складання. Це порівняно просто виконати у залізі (розписувати вже не буду).

(Чи не вважати FMA «простою» операцією? Та хай, воно ж не сінус...)

А ще є підходи, які дозволяють підвищувати точність операцій через декілька змінних і корекцію в старших змінних; для додавання і віднімання прикладом є Kahan summation, але за аналоґією можна це робити також для інших операцій (значно складніше).

2. Можуть залишатись засоби, які рахують по-старому або інакше.

В IBM S/360 і проміжне, і фінальне округлення робилось просто усіченням (округленням до нуля), тримали тільки одну сторожову цифру (guard digit), у цій системі шістнадцяткову, для того, щоб при зсуві вліво (як показано в частині 3 для віднімання) зовсім точність не спаплюжилась. Більш того, в початкових версіях в S/360 ця guard digit не була в реалізації з подвійною точністю, бо автори вирішили, що сама по собі подвійна точність вже достатня; це настільки псувало результати (не самою по собі точністю, а її нестійкістю, на додачу до нестійкости, що викликана основою 16), що в 1968-му всі екземпляри S/360 були перероблені (microcode update, в сучасних термінах) — додали цю цифру. Усічення залишилось, але вже мали дещо стабільнішу поведінку.

Зараз всі актуальні компʼютери лініі SystemZ підтримують і IEEE754 (не тільки двійкову версію, а і десяткову). Але кожен, хто з ними працює, має ненульовий шанс зустріти таке леґасі (зветься «HFP»).

В PDP-11 FPP проміжнє округлення було усіченням, хвіст був 1 біт для одинарної і 3 для подвійної точности, перевірявся тільки старший біт хвоста; фактично це виконувало щось віддаленно схоже на half-away (або toward zero при відповідному біті в режимі сопроцесора), але без повної точности. Це вже середина 1970-х, і все одно розробки були буквально на рівні камʼяного віку... десь тоді всім и припікло.

Про X86 FPU (X87) вже було сказано: з точністю у нього гірше, ніж у проґрамних реалізацій, і, що головне, неконтрольовано. (А проґрамна реалізація будь-чого складніше чотирьох операцій арифметики там буде значно дорожча, ніж з SSE, через застарілий метод стека регистрів.) І це стосується не тільки складних функцій, як sin(); окрему проблему вносить його внутрішній формат даних. Про це буде окрема частина. Скоріш за все, читачам цього вже не доведеться стикнутись з ним: майже всі сучасні засоби використовують в 64-бітному режимі тільки SSE (або інші архітектури), а 32-бітний X86 залишився тільки в специфічних використаннях.

Але задля історії... ось подивімось на деякі деталі реалізації в Pentium. Не тільки внутрішній формат з 64 бітами мантиси, а ще й 67-68 біт на обчислення одної команди, і навмисно спотворені коефіцієнти ряду Тейлора для зменшення похибки у реальних умовах. (Любителям «технопорну» — «саме воно», і ще один блоґ.)

З x87 FPU може бути ще одна проблема: проміжне представлення чисел в «розширеному» (extended) форматі дає інші результати, ніж одразу в цільових (як double). Про це буде окреме згадування.

Округлення

Головна тонкість вже згадувалась в частині 3: мало які платформи реально дозволяють контролювати режими округлення. «Слухай, Петрику, пісню про комбайн», тобто використовуйте half-even і не думайте про альтернативи.

Відсутність RPSP в стандарті дивує тим, що, з одного боку, приховується фундаментальний режим, який і так реалізований в залізі, а з іншого, не дає акуратно контролювати тонкі випадки, коли, наприклад, виконуємо операції в double з фінальним скороченням точности до single, а це не рідкий варіант, наскільки я бачу з прикладів коду. Проте для звичайної двійкової арифметики (перший домен зі вводної частини) це нечасто створює проблеми, а там, де вони є, краще одразу виконувати все в double (чи яка там буде потрібна точність) і не втрачати на конверсії... так що тут економія «на сірниках», здається, не створює надвеликих проблем.

Про леґасі не-IEEE754 варіантів дивіться вище.

З прикрих курйозів треба згадати round() в Java (до цілого): конкретно, java.math.Round(). Зазначений для нього режим округлення — half-ceil (half-to-plus-infinity); він не має жодного реального застосування, чи я його ніколи не бачив. Обґрунтування такого вибору я не знайшов. Тому зробити в Java округлення до цілого за найбільш типовим half-to-even вимагає трюків. Стара round() із C визначає half-away, з нею теж не збігається.

Так мало того, Javascript повторив це half-ceil! Інакше як гонитвою за імітацією, цю маячну не можна пояснити.

Округлення на краю нескінченности, направлені режими

Один з досить неочікуваних ефектів виглядає так: візьмемо операцію, яка має привести до переповнення — наприклад, 1e308 * 1e308 в double. При режимах округлення до нуля (toward zero) або донизу (toward -∞), результатом буде не переповнення з видачею Infinity, як при дефолтном half-to-even, а максимальне скінченне з double (1.797...e308). На це є явна вимога стандарту. Стандарт каже (пункт 4.3.2): «With roundTowardNegative, the result shall be the format’s floating-point number (possibly —∞) closest to and no greater than the infinitely precise result». Дзеркально таке ж саме спостерігається біля мінус нескінченности.

Можна зрозуміти цей підхід і результат, але проміжок між таким максимальним числом і нескінченністю дуже великий 😱

Два нулі

В кожному з форматів можна представити +0.0 і —0.0, це різні значення. При цьому порівняння скаже, що вони рівні одне одному... але вони не однакові за результатом деяких операцій, наприклад:

>>> math.atan2(+0.0, +0.0)
0.0
>>> math.atan2(-0.0, +0.0)
-0.0
>>> math.atan2(+0.0, -0.0)
3.141592653589793
>>> math.atan2(-0.0, -0.0)
-3.141592653589793

atan2 вираховує кут полярних координат точки з показаними декартовими координатами (в порядку: y, x). Тобто, якщо ми показуємо запит на визначення направлення вектору, який дивиться з початку координат у напрямку вздовж x наліво, то кут дорівнює π, а якщо направо — нуль. Але і π тут два різних. (Звісно, в визначенні є деяке свавілля стандартизаторів; два нулі могли б показувати будь-який кут у своєму квадранті, хоч 45, хоч 23.4567 градусів. Але і воно має свої причини.)

Отримати —0.0 можна або вирахувавши новий «чесний» 0 при округленні до -∞, або в діях типу скласти два —0.0. Не залізатиму тут в деталі. Для чого він взагалі був потрібний?

Тут склались різні потреби. Одна — якщо замале значення втрачене повністю і ушло в 0, то хоч знак буде видний. Другий — тонкі ефекти особливих функцій, як вище про atan2().

Звісно, це теж викликає обурення і голоси, що так не треба було робити. Але обміркування показує, що по більшости таки від двох нулів є користь. Проблема виникає на границях підходів. Якщо скласти два значення, +0.0 і —0.0, сумма буде —0.0 тільки при округленні до -∞. При інших, включаючи найзвичніший half-even, вона дорівнюватиме +0.0. Чи коректно це? А у нас є вибір? Нітъ, бо третього варіанту нуля, без знаку, немає. Ну, хоч для більшости випадків воно вирішено.

Звісно, це далеко не все, що можна було б назвати проблемами стандарту. Далі буде, наприклад, про експорт-імпорт і порівняння.

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