Rust std и musl неверно округляют fmaf на редких числах

Разработчик под ником Shnatsel 2 сентября описал, как, реализуя fused multiply-add для библиотеки fearless_simd, нашёл ошибку округления на субнормальных числах, а затем ту же ошибку обнаружил в f32::mul_add стандартной библиотеки Rust, в std::simd и в функции fmaf() из libc musl. Результат отличается от правильного на один младший бит; в исходном контрпримере это около 15 частей на миллион.

Ошибка проявляется только там, где FMA эмулируется программно: на x86 это дешёвые или старые Intel без AVX2 и процессоры Hygon без аппаратного FMA, а также 32-битный ARM и старые встраиваемые тулчейны. 64-битный ARM, по словам автора, не затронут. На практике это касается детерминированных симуляций и численных библиотек, где одинаковый код обязан давать одинаковые биты на разных машинах: сетевые игры с lockstep, физические движки, воспроизводимые вычисления. Автор честно пишет, что не знает, насколько ошибка важна в реальных программах.

Что такое FMA и откуда берётся лишний бит

Fused multiply-add вычисляет a*b+c с одним округлением в конце вместо двух: после умножения и после сложения. Так требует IEEE 754, и именно на это рассчитывают численные алгоритмы: одно округление даёт предсказуемую ошибку. Когда у процессора нет аппаратной инструкции, библиотека эмулирует FMA через операции с расширенной точностью и ручную коррекцию округления. Ошибка, которую нашёл автор, живёт в этой коррекции на субнормальных числах, то есть очень маленьких значениях, у которых часть мантиссы уже ушла в ноль. В таких случаях программная реализация округляла не в ту сторону.

Автор переносил формально верифицированный алгоритм и гонял его в CI на эмуляторах SIMD, добавив тесты на миллион субнормальных значений и ещё миллион случаев, требующих коррекции округления. Отдельно он попросил языковую модель искать контрпримеры, и она нашла случай, ломающий первую версию исправления для стандартной библиотеки Rust: там ошибка приводила ещё и к неправильной установке флагов исключений плавающей точки.

Как проверить свою платформу

Достаточно трёх чисел. Ниже проверка на C для fmaf из системной libc; в Rust то же самое делается через f32::from_bits и mul_add:

			#include <math.h>
#include <stdio.h>
#include <string.h>

static float bits(unsigned u) { float f; memcpy(&f, &u, 4); return f; }

int main(void) {
    float r = fmaf(bits(0x97000800u), bits(0x1cfff001u), bits(0x00010002u));
    unsigned out; memcpy(&out, &r, 4);
    printf("%08x (ожидается 00010001)\n", out);
    return 0;
}
		

Если вывод 00010002, ваша libc или тулчейн эмулируют FMA с ошибкой. На машине с аппаратным FMA (любой современный x86 с AVX2, 64-битный ARM) результат будет верным независимо от библиотеки, поэтому проверять имеет смысл именно сборки под старое или встраиваемое железо и контейнеры на musl (Alpine) на таких хостах.

Состояние исправлений

• fearless_simd: исправлено, можно использовать как реализацию FMA без этой ошибки.
• Rust, тип f128 на 32-битном ARM: ошибка ABI исправлена в компиляторе (тип доступен только в nightly).
• Стандартная библиотека Rust (f32::mul_add, std::simd): патч на момент публикации ожидал ревью; номер версии с исправлением не назван.
• musl: предложены минимальное исправление и полная переработка fmaf(), ни то ни другое не влито. Автор отмечает, что код FreeBSD, откуда растёт реализация, мог копироваться и в другие тулчейны.

Кому это ничего не меняет: сервисам на современных x86 и ARM64, где FMA аппаратный, и коду, не зависящему от битовой точности. Кому стоит проверить: авторам физических движков и симуляций с детерминизмом между машинами, разработчикам под 32-битный ARM и тем, кто собирает статические бинарники под musl для старых серверов. Про то, как выбор libc влияет на скорость тех же бинарников, мы писали на примере Rust Coreutils; смежная тема из этой же недели: переписывание линкера mold на Rust. Что мы проверим дальше: когда патч попадёт в стабильный Rust и появится ли исправление в релизе musl.

Источники: Shnatsel: Implementing FMA and finding bugs in C and Rust standard libraries, Обсуждение на Lobsters