Плаваюча (рухома) крапка, частина 8: різне. Субнормальні. Розширений 80-бітний формат. Баґ 323

Частини 1-2: ввідна, перші граблі.

Попередня: частина 7: десяткова рухома крапка.

Ця частина мала і стосується ефектів, які навіть на тлі тем цього циклу дуже спеціалізовані :)

Субнормальні (денормалізовані) числа

Почнемо з такого прикладу. В частині 7 були розглянути числа з основою 10 і показано, що для них примусова нормалізація не дає суттєвої переваги (хоча це не 100% твердження; для когось відсутність необхідности зсуву більшого числа, як показано в виконанні операцій, може мати значний сенс). Для таких чисел когорти представлень містять декілька рівнозначних варіантів, як, наприклад, при 6 цифрах 123000e3, 012300e4, 001230e5, 000123e6. У реальному представленні порядок завжди має обмеження. Нехай він може бути від −99 до +99. Якщо у нас число представлене як 000123e-96, воно ж представляється 123000e-99. Але для 000012e-96 вже не буде представлення зі старшою цифрою не 0 (тобто нормалізоване); найкраще, що ми тут отримуємо, це 012000e-99. Це і є «субнормальне» (subnormal) число.

В сучасній термінолоґії кажуть саме subnormal (субнормальне), але багато де застиг термін denormalized (денормалізоване) або іменник denormal (наприклад, SSE флаг «flush to zero»).

В двійкових форматах IEEE754 в основному варіанті (див. частину 3) є прихована старша 1 у мантиси. Це зручно, якщо нормалізація форсована; якщо ми маємо завжди старшою цифрою 1, можна її не зберігати і тим самим виграти аж 1 біт. Звідси витікає очевидний підхід: якщо ми універсалізуємо таке представлення, субнормальних вже не буває. Реалізації рухомої крапки до IEEE754, що вимагали нормалізації представлень, зазвичай такі числа перетворювали в нуль. Але одна з найцікавіших рис IEEE754 в тому, що він намагається зберегти такі значення якомога точніше, за рахунок невеличкої додаткової складности реалізації.

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

Наприклад, у 16-бітному форматі (5 біт на порядок і 10 на мантису не включаючи прихований біт):

0:00001:0000000000 — найменше нормалізоване число, дорівнює 2↑-14 == 1/16384 ≈ 6.1035e-5;
0:00000:1111111111 — найближче до нього (найбільше) субонормальне == (2↑-14 - 2↑-24) ≈ 6.0976e-5;
… продовжуємо вниз …
0:00000:0000000001 — найменше за модулем число, що може бути взагалі представлене і не нуль = 2↑-24 ≈ 5.96e-8;
0:00000:0000000000 — +0.0 (представлення з порядком субнормального, особливий випадок).

(Можна переінтерпретувати наймолодший біт поля порядку як цей найстарший біт мантиси, якщо решта бітів порядку дорівнює 0.)

При ґенерації таких чисел має місце поступова (gradual) втрата точности (зріст відносної похибки): для подвійної точности відносна похибка буде 1/2↑53, 1/2↑52, ..., 1/4, 1/2 для найменших значень. Якщо результат обчислень став субнормальним і при цьому округленим, стандарт вимагає підняття одночасно двох флаґів помилок: Underflow і Inexact; якщо число отримане без округлення, то флаґи помилок не змінюються. Контролюючи значення Underflow, можна перевіряти на таку поступову втрату точности. Звісно, якщо вам такий флаґ дали (знаємо, що стандартні реалізації майже всюди, крім C/C++, не дають доступу до такого стану).

Що не так? А ось тут є цікавий і сумний момент. Якщо подивитись у реалізацію Berkeley Softfloat (згадував вже її), підтримка субнормальних, звісно, чогось коштує, але це максимум додатково 1-2 прості операції на одне обчислення. Ну да, треба це зробити дуже акуратно; коли отримуєш проміжний результат і його треба зсунути вправо так, щоб умістилось як раз при субнормальному зберіганні і дало потрібне округлення... навряд чи хтось з першого разу зробить це правильно, і потрібні тести для підтвердження вичитки коду. Але у апаратних реалізаціях IEEE754 до якогось моменту, у Intel x86 (FPU і початкове SSE), Motorola... субнормальні приводили до виходу мікропроґрами на «повільний шлях» (slow path), де на одну операцію могли додаватись сотні тактів! Я не знаю, як і чому вони написали ці мікропроґрами так складно. В Internet нема (поки що?) опису цих блоків з деталями реалізації, максимально, що є — це «ну вибачайте, випадок тут особливий, нетиповий і тому неоптимізований».

Через ці дивнощі в x86 SSE були введені флаґи режиму: flush-to-zero (FTZ) і denormals-are-zeros (DAZ) в реґістрі MXCSR, які просто вимикали ці довгі гілки виконання. Для ARM є біти в FPCR: FZ, AH, NEP; AArch32 в векторних командах безумовно перетворює субнормальні на нулі (хм).

Що ж зараз? Як кожен може помітити, після, приблизно, 2005 року (для x86 це кінець лінії Pentium 4, успіх Core), зупинився рост тактової частоти і почалось удосконалення вже глибокої лоґіки виконання операцій. Перелом пройшов, для Intel, на поколінні Nehalem. І Intel з AMD одночасно почали змінювати реалізації... але різними шляхами. Вже в поколінні Sandy Bridge в SSE (не FPU) Intel усунув значну частину непередбачуваних затримок з субнормальними; приблизно тоді ж AMD, навпаки, почав усувати їх з FPU (не SSE). Зараз у актуальних моделях обох головних виробників x86, наскільки мені відомо, таких затримок немає. (Я не знаю деталей апаратної реалізації, які до того ж у кожної моделі можуть бути свої, тому не знаю, як все перевірити. Але тести на базову аріфметику показують, що в основних операціях проблема усунута. Можете перевірити самі.)

Чи означає це, що підіймати FTZ і DAZ (чи їх аналоги у інших архітектур) зараз не треба? Мабуть, не завжди означає: у щільних шляхах виконання кожний такт може мати значення, і якщо ви впевнені, що такі числа не мають сенсу у ваших операціях, то можна і підняти їх. Але якщо такого нема, то краще залишити неактивними і отримати більше точности.

Ще один супер-тонкий момент з субнормальними це коли вони визнаються такими у погранічних випадках. З тим же десятковим прикладом з початку частини, якщо у вас до округлення число дорівнює 99999.998e-99, а округлення дасть зміну на більше (100000e-99, яке вже нормалізоване), то буде число детектоване як субнормальне чи ні? І тут стандарт додає ще специфіки:

The underflow exception shall be signaled when a tiny non-zero result is detected. For binary formats, this shall be either:
a) after rounding — when a non-zero result computed as though the exponent range were unbounded would lie strictly between ±b↑emin, or
b) before rounding — when a non-zero result computed as though both the exponent range and the precision were unbounded would lie strictly between ±b↑emin.
The implementer shall choose how tininess is detected, but shall detect tininess in the same way for all operations in radix two, including conversion operations under a binary rounding attribute.
For decimal formats, tininess is detected before rounding — when a non-zero result computed as though both the exponent range and the precision were unbounded would lie strictly between ±b↑emin.

Примітимо, що сам по собі факт округлення ґенерує Inexact; а ось чи видасть така реалізація флаг Underflow, якщо вихідний результат не субнормальне, а нормальне (це може бути мінімальне з нормальних) — залежить від реалізації. Більшість таких для двійкової арифметики використовує after rounding, бо інакше встановлення флаґа Underflow одночасно з нормальним результатом надто бентежитиме; але фіг ви знайдете це в документації на процесор, це тільки шукати інші джерела або самому експериментувати. А ось для десяткової детект before rounding вже встановлений стандартом, бо так точніше(?)

У любому випадку, це настільки «тонка» матерія, що я активно сумнівався, чи згадувати це взагалі. Фінально вважаю, що приклад «з чим доводиться стикатись у погранічних випадках» тут вартий пів-сторінки.

«Розширений» формат точности x87 і Motorola

Те, що описане в цьому і наступному підрозділі, щодалі отримує суто історичне значення, але, якщо хтось використовує 32-бітний режим x86 (ще може бути) або явно просить FPU, то для нього це має реальне втілення. Описуватиму x87 (тобто, 8087, 80287 і далі, поки не стало зватись просто FPU в середині процесора); для Motorola деталі інші, але головна ідея збігається. (Ця лінія Motorola вимерла в процесорах загального використання, але є в деяких embedded рішеннях... ще є шанс зустріти.)

Доступність FPU залежить від компілятора. Звісно, на асемблері можна написати що завгодно. В більшості сучасних юніксів для x86, в 64-бітному режимі замовчування — SSE; щоб використати FPU, треба прийняти особливих заходів (для GCC, Clang це опція -mfpmath=387). З Microsoft C/C++ це неможливо взагалі в 64-бітному режимі.

Внутрішнє представлення в FPU завжди 80-бітне, з яких 1 — знак, 15 — порядок і 64 — мантиса (significand), яка явно включає в себе старший біт (на відміну від стандартних транспортних форматів IEEE, де він прихований). Читати з памʼяти і писати можна в одному з трьох варіантів (32, 64, 80 біт). Обчислення в FPU завжди виконуються в діапазоні порядків як для 80-бітного режиму (тобто 15 біт на зсунутий порядок => діапазон від −16383 до +16384, субнормальні як звичайно), але точність мантиси можна обмежити режимом (24, 53, чи 64 біта — відповідно, для 32-, 64- чи 80-бітного формата).

(Чому всередині використовується такий розширений формат? Історія каже, що це був єдиний доступний на кінець 1970-х метод зробити ґарантовану точність на складних трансцедентних операціях, ба навіть на «простому» (хм) піднесенню в степінь, для якої не було одної команди, а був хитрий ланцюжок: фактично, замість x↑y вираховувалось 2↑(y*log2(x)), при чому 2↑x виконувалось через більш точну 2↑x-1, а log2(x) — через log2(x+1)...

В деяких випадках можна це вміння FPU використовувати для ще більшої точности, ніж звичайне double. Але реалізація деяких функцій вимагатиме окремих бібліотек, бо у стандартних не буде адекватної точности на всі 64 біти.

Вже потім IEEE754 стандартизував розширені (extended) формати для кожного з базових форматів, з тих же причин: підтримка підвищеної точности. Але я не знаю апаратну реалізацію після 2000-го, яка б їх розвивала. Це вже викликано переносом забезпечення точности складних операцій на софт.)

Де 80-бітний формат доступний, він зветься в C/C++ як long double. Де він недоступний, це long double є аліасом для просто double. Ще додає специфіки те, що в типовій компіляції numpy (Python) тип, який зветься в ньому float128, насправді реалізує цей 80-бітний формат; але і на це не можна строго закластись, треба перевіряти конкретну реалізацію.

Я не вдаватимусь у деталі реалізації FPU, таких, як стекова модель реґістрового пулу, взаємодія з основним процесором... — все це за межами нашої теми. Але є дуже суттєвий наслідок такої реалізації: «bug 323».

«Баґ 323»

Це дивне явище має місце саме з FPU і через його розширений формат. Названий номер це з GCC bugzilla, але став універсальною назвою проблеми. (Аналоґічно, «bug 12309» всесвітньо відомий як безглуздий підхід до керування памʼяттю в Mach-like VM, особливо і особисто в Linux.) За посиланням можна знайти приклад дивної поведінки проґрами:

void test(double x, double y)
{
  const double y2 = x + 1.0;
  if (y != y2) printf("error\n");
}

void main()
{
  const double x = .012;
  const double y = x + 1.0;

  test(x, y); // може видати "error"
}

Або приклад того ж в C#, ще ясніше — передаю як було в джерелі:

    static void Main(string[] args)
    {
        // !! Саме так, дві однакові операції!
        float f = Sum(0.1f, 0.2f);
        float g = Sum(0.1f, 0.2f);
        // Ефект залежить від закоментованости наступного рядка...
        //System.Console.WriteLine("f = " + f + " and g = " + g);

        if (f == g) { System.Console.WriteLine("The sums are equal"); }
        else { System.Console.WriteLine("The sums are Not equal"); }
    }
    static float Sum(float a, float b)
    {
        System.Console.WriteLine(a + b);
        return a + b;
    }

В Debug в тестах, де це було знайдено, видає equal, а в Release — Not equal.

В цій формі навіть «інтерн» побачить, що щось не те. 😲😠

В чому проблема? Якщо значення f і g одночасно збережені в реґістрах FPU, результати будуть однакові. Якщо обидва записані в памʼять, те ж саме. Але якщо тільки один з них вивантажений в памʼять, а потім назад завантажений в FPU, то це вивантаження в памʼять приводить до округлення з розширеного (64 біта мантиси) до одинарної точности (24 біта), яке змінює значення; і тоді порівняння дає нерівність. Але якщо в режимі роботи FPU була встановлена одинарна точність для кожної операції, ефект, звісно, зникає. В Windows за замовчуванням налаштована в FPU подвійна точність; в Linux — розширена. Для реґулювання є функції з іменами типу _control87(), _controlfp(), fpsetprec().

Ще один варіант ефекту.

Для проґрам на C має бути визначений FLT_EVAL_METHOD, але, якщо він дорівнює −1, це значить «нічого не ґарантуємо» :( В Linux/x86 в 64-бітному режимі він дорівнює 0 (нема розширень типу і відповідно нема таких ефектів), а в 32-бітному — 0 при SSE, і 2 при FPU (все через long double, ефект в повний зріст). GCC, Clang знають флаґ -ffloat-store, який вимагає збереження всіх результатів операцій, включаючи проміжні, в памʼять; це фіксить проблему ціною виробности.

Основна частина ефекту «баґа 323» — емоційна. Якщо пишеш код, очікуєш стабільности у деяких базових характеристиках, як ідемпотентність операцій і незмінність отриманих значень. І саме тут при зіткненні з даним ефектом виникає відчуття повної відсутности стабільного ґрунту під ногами. Може, це додало свій внесок в уход від FPU...

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