Введённый в 1985 году стандарт IEEE-754 для чисел с плавающей запятой был предназначен для решения проблемы разнородности реализаций чисел с плавающей запятой, мешавших портируемости кода, а также для повышения стабильности между платформами.

Он получил широкое применение и многократно пересматривался в течение прошедших лет. Если вы когда-нибудь работали с любыми вещественными числами в своих приложениях, то они, вероятно, отвечали этому стандарту.

Моя работа в течение последнего года заключалась в анализе погрешности различных математических функций, накопления этой погрешности и способов её уменьшения при помощи различных программных паттернов. Одной из исследованных мной тем были базовые математические функции, используемые в функциях активации нейронных сетей, а также способы их аппроксимации для повышения производительности. В процессе работы нам пришлось столкнуться с противодействием со стороны людей, активно стремящихся к корректной реализации математических функций и к соответствию их стандартам, в частности, к соблюдению обеспечения корректности одной наименее значимой единицы измерения (unit in last place, ULP) для элементарных функций.

Я был заинтересован в дальнейшей работе по аппроксимации этих функций, поэтому приступил к исследованию того, каким образом они гарантируют корректность, и если они корректны только на 1 ULP, то где располагаются ошибки в области определения функции.

В процессе изучения я обнаружил, что ни одна из популярных математических библиотек, используемых во множестве сфер вычислений, на самом деле не выполняет корректное округление в соответствии с требованиями любой версии IEEE 754 после первой редакции 1985 года.

▍ Реализация


В качестве примера я возьму функцию sin, потому что её относительно легко объяснить, и она хорошо продемонстрирует все аспекты, которые я стремлюсь показать.

Функция sin — это периодическая симметричная со смещением функция относительно оси Y, каждый сегмент кривой имеет при каждом повороте одну и ту же форму. Кривая синуса повторяется каждые 2pi, а часть кривой с отрицательным Y является отражением части кривой с положительным Y. И каждую кривую в кривой синуса можно построить из её части.

Это позволяет нам сделать две вещи: воспринимать любой x во вводе как остаток от x/2pi, чтобы понять, в какой части интервала 0-2pi мы находимся (это называется снижением размерности), а затем определить, в какой части кривой мы находимся, и обрабатывать более точную аппроксимацию этой более простой кривой.

Если не вдаваться в подробности, то реализация выглядит примерно так:

function sinnaive(x::AbstractFloat)
    # Получаем значение и снижаем размерность до [0,pi/2]
    C     = pi/2.0
    k     = Int64( floor( x /C ) )
    phase = floor(k) % 4
    
    # наивное снижение размерности
    y = x - (k * C)
    
    if(phase == 0)
        return  approxsincurve(y)
    elseif (phase == 1)
        return  approxcoscurve(y)
    elseif (phase == 2)
        return -approxsincurve(y)
    elseif (phase == 3)
        return -approxcoscurve(y)
    else
        @assert false ("Invalid Range Reduction")
    end
    
    return Inf
end

В этом и состоит секрет получения очень хороших результатов синуса, зависящих от получения очень хорошей аппроксимации частей кривой и проверки того, что мы находимся в нужном месте. Существует ещё множество разных трюков, делающих подобные реализации идеальными. В основном они связаны с обеспечением корректного округления подвыражений и использованием по возможности команд с высокой внутренней точностью (например, Fused Multiple Add).

Но поскольку всё это на базовом уровне зависит от аппроксимации, очень сложно обеспечить возможность получения идеально точного ответа при работе с ограниченной точностью. И из этого мы получаем предел погрешности в 1 ULP, о котором говорят все инженеры, работающие с такими функциями.

Но здесь есть любопытный момент: так как конкретная техника, которую нужно использовать, в стандарте не определена, а корректность может зависеть от аппаратных команд, в результате у нас есть различные реализации, дающие погрешность в разных позициях. Например, вот функция sin из библиотеки Microsoft CRT, которая применяется для C и C++ в компиляторе MSVC, и результаты реализации sin на языке Julia (позиционируемом как быстрая альтернатива Python в числовых вычислениях).



Каждая из вертикальных линий обозначает некорректно округлённое значение. Суммарно у Julia их меньше, чем в реализации Microsoft. Однако они гораздо сильнее разнесены; Microsoft же выбрала реализацию, приводящую к высоким уровням корректности в больших блоках, а все погрешности сконцентрированы в областях примерно через каждые pi/2.

Если придерживаться того, что говорили нам инженеры, то обе реализации точны в пределах 1ULP, но добиваются этого очень по-разному. При таком способе реализации мы, по сути, создаём сторонний канал информации, позволяющий пользователю точно узнать, с какой библиотекой выполнялась компоновка, если он может запрашивать в вашем приложении базовые функции. Вероятно, по большей мере эта информация безвредна, но обычно о её наличии большинство пользователей не догадывается.

▍ Несоответствие стандарту


Именно мысль о стороннем канале утечки информации заставила меня обратиться к различным сообществам разработчиков математического ПО, чтобы понять, знают ли люди об этом факте, и не приводили ли различия в концентрациях погрешностей к каким-то проблемам. Я рассуждал так: «Если погрешность сконцентрирована, то, вероятно, это приводит к меньшей согласованности, и алгоритмы градиентного спуска могут иметь пониженную эффективность, поскольку последовательные числа с большей, чем средняя, вероятностью будут одинаковыми при округлении не в том направлении».

Меня уверили, что абсолютная погрешность в числах слишком мала, чтобы стать проблемой. Честно говоря, я этому не поверил. Я работал над играми, в которых даже крошечные погрешности вычислений с плавающей запятой превращались в визуальный хаос. Поэтому я обратился к стандарту, чтобы проверить, есть ли в нём какие-то дополнительные требования к этой погрешности в 1ULP, чтобы избежать различий между платформами, для чего изначально и был создан этот стандарт.

И тут я выяснил, что ограничение в 1ULP для элементарных функций в стандарте отсутствует. Его нет в версии 1985 года: она не накладывает никаких ограничений на корректность этих функций. Насколько я понимаю, единственное их упоминание заключается в том, что их разработку следует вести на основании стандарта.

В следующей ревизии IEEE 754-2008 в этом отношении многое изменилось. В ней есть «Section 9. Recommended operations»; в этом разделе приведён список всех рекомендуемых функций, исключения, которые они должны выбрасывать, их область определения и, что самое важное — требуемая степень их корректности:

«Отвечающая стандарту функция должна возвращать результаты, корректно округлённые в применяемом направлении округления, для всех операндов в области её определения».

Что же означает «корректно округлённые»? Это означает погрешность менее чем в 0,5 ULP. То есть должен быть ноль несоответствий между корректным результатом и результатом, возвращаемым функцией. Иными словами, все версии математической библиотеки, отвечающей стандарту IEEE, должны возвращать одинаковые результаты в одинаковом режиме округления. Это утверждение сохранилось и в IEEE 754-2019, выпущенном в 2020 году.

Все протестированные мной стандартные математические библиотеки, по утверждению разработчиков соответствующие стандарту IEEE, не удовлетворяют этому ограничению.

▍ Реакция


Сначала я поговорил со старым приятелем, с которым мы когда-то работали в компании-разработчике игр, потому что у нас возникали проблемы с этими функциями, когда мы пытались реализовать детерминированную генерацию планет. У пользователей на разных машинах генерировались разные паттерны, то есть это ломало многопользовательский режим нашей игры. Выслеживая ошибку, мы добрались до GPU, имеющих собственные правила, и до погрешности в 1ULP на стороне CPU. Поэтому я и захотел начать раздел с этой истории из жизни, показывающей, что в некоторых приложениях такое несоответствие стандарту имеет реальные последствия.

Затем я поделился своими открытиями с сообществом разработчиков на Julia (потому что я использовал этот язык для численных вычислений в своих инструментах анализа аппроксимаций) и с сообществом C++, потому что с этим языком у меня больше опыта.

Реакция на мою информацию оказалась очень оживлённой, но обескураживающей. Похоже, многие люди, участвовавшие в реализации этих функций, осознают, что они не отвечают требованию стандарта, но закрывают на это глаза. Другие делали заявления, доказывающие, что они не понимают или не читали стандарт. Ещё кое-кто полностью отвергал требования стандартов, мотивируя это их чрезмерной строгостью.

Приведу очень интересные цитаты из ответов на мой вопрос о том, почему об этом знают, но не исправляют. Кто-то признавал, что не соблюдает стандарт, другие возлагали бремя обеспечения корректности на пользователя, который должен знать, что функции ошибочны:

«Функции соответствуют спецификации, которой, по мнению людей, они соответствуют».

«Думаю, что есть более масштабная проблема: многое из этого не заявляется явным образом; к тому же рекомендации спецификации IEEE чудовищно непрактичны».

«Если я пишу код, который поломается, если вызываемые им функции обеспечивают гарантию точности только в 1 ULP, то такой код сам по себе похож на карточный домик, и достаточно лёгкого дуновения, чтобы он развалился. Мне следует просто использовать числа с повышенной точностью, а не требовать от мейнтейнеров библиотеки рушить производительность у всех пользователей, косвенным образом используя высокую точность даже в тех ситуациях, когда они не понимают, важно ли это».

«В конечном итоге, всё сводится к затратам: стоят ли эти лишние 0.000897 ULP точности затрат производительности, которых они требуют? (а также усилий для корректной реализации, то есть для решения достаточно нетривиальной задачи)».

Однако самый понравившийся мне ответ я получил от сообщества Julia. Он напоминает вызов:

«Автор, а ты сам собираешься отвечать на все жалобы пользователей, желающих знать, почему Julia в тридцать раз медленнее, чем все остальные числовые системы?»

Я не такой уж хороший писатель и не самый дипломатичный человек. Но меня волнует производительность, и мне кажется, что тридцатикратное замедление при корректной реализации округления 32-битных чисел — это преувеличение. Поэтому я решил, что решу проблему жалоб, реализовав округление так, чтобы оно не было медленнее в тридцать раз. Жалоб не будет, а значит, моя работа выполнена!

▍ Корректная реализация согласно IEEE


Оказалось, что у большинства современных 64-битных процессоров есть достаточно мощные модули для выполнения операций с 64-битными числами с плавающей запятой. Поэтому я подумал: почему бы нам просто не реализовать 32-битный sin как 64-битный sin с корректным округлением? Теоретически в большинстве случаев округлённое значение будет корректным; некорректным оно будет в случае плохого двойного округления, которое может возникать, но при плотности чисел в интервале вывода от -1.0f до 1.0f это и чрезвычайно маловероятно, и крайне легко проверить, потому что на современной машине мы тривиальным образом можем протестировать каждое входное значение в интервале сниженной размерности от 0.f до 2.f*pi и даже пойти дальше, если нам нужна абсолютная уверенность.

Итак, мы реализуем 32-битный sin следующим образом:

function IEEE_sin(x::Float32, correctlyRounded=true)::Float32
    if(correctlyRounded)
        return Float32.(sin(Float64(x)))
    else
        return sin(x) 
    end
end

function IEEE_sin(x::Any, correctlyRounded=true)::Any
    return sin(x) 
end

А затем протестируем его.

Актуальная реализация в Julia:

sin: Max Err(ULP): 5.00897385615018681127840279120244109e-01. Mean Err(ULP): 2.50005459168491000625254536923056698e-01
sin: Max Err(Abs): 5.9604645e-8. Mean Err(Abs): 7.545918e-12
sin: Количество несоответствий: 10827

Наша реализация, соответствующая стандарту IEEE:

IEEE_sin: Max Err(ULP): 4.99999995510513399926379631449132211e-01. Mean Err(ULP): 2.50005375827442751842250455717785186e-01
IEEE_sin: Max Err(Abs): 0.0. Mean Err(Abs): 0.0
IEEE_sin: Количество несоответствий: 0

Сравнение производительности:

Base.sin: Trial(1.042 ms)
IEEE compliant sin: Trial(1.041 ms)

В результате мы получили функцию, выполняющую корректное округление, согласованную в разных реализациях и прогонах (судя по бенчмаркингу с помощью стандартных инструментов Julia) с примерно теми же затратами.

Я повторил эти тесты с большой выборкой рекомендованных IEEE 754 функций. Этот подход не сработал для log, но количество несоответствий было близко к нулю. Существуют разные способы исправить это. Медленный способ — это повышение точности, но я также экспериментирую и с измерением расстояния до ближайшего float при повышенной точности с последующим выбором. Похоже, это решает проблему.

Под конец проделанной мной работы у меня была выборка из большинства рекомендованных функций с корректным округлением в 32 битах. Это сработало на всех протестированных мной языках.

▍ Заключение


Подведём итог: похоже, мы можем обеспечить соответствие стандарту IEEE 754 в 32-битных реализациях на современных машинах без существенного ущерба производительности. Корректность для всех значений при повышенной точности проверить более затратно и сложно (table maker's dilemma), но это не причина закрывать доступ к этой согласованной со стандартом функции вместо того, чтобы выбирать произвольный уровень погрешности и придерживаться его без ограничений отрицательных побочных эффектов, к которым это может привести.

В конечном итоге, смысл стандарта — предотвратить неожиданное поведение, а в современных реализациях есть большая доля неожиданного поведения. Можно начать с исправления 32-битных реализаций, а затем перейти к более сложным задачам.

Меня критиковали за то, что я якобы не хочу, чтобы существовали быстрые решения с 1ULP. Я не против них: на самом деле, я хочу, чтобы динамическая библиотека реализовала произвольные уровни точности и мы не тратили вычислительные ресурсы на ненужную работу. Похоже, все мы неосознанно делали это, потому что все библиотеки втайне решили не соблюдать эту часть стандарта.

Однако в своей динамической библиотеке я обеспечу детерминированность результатов между реализациями. sin(x) в библиотеке А должен совпадать с sin(x) в библиотеке Б. Без фундаментальных идей эквивалентности функции и значений вся модель вычислений с плавающей запятой разваливается непредсказуемым и пугающим образом.

Telegram-канал со скидками, розыгрышами призов и новостями IT ?

Комментарии (36)


  1. SuAlUr
    12.02.2025 13:23

    Фундаментально!


  1. unreal_undead2
    12.02.2025 13:23

    похоже, мы можем обеспечить соответствие стандарту IEEE 754 в 32-битных реализациях на современных машинах без существенного ущерба производительности.

    Похоже, тестировался невекторизованный код.


  1. kibb
    12.02.2025 13:23

    Статья 5ти летней давности. Вот что-то не такое старое https://members.loria.fr/PZimmermann/papers/accuracy.pdf


    1. Balling
      12.02.2025 13:23

      Вы главное забыли упомянуть, в этой статье упомянуто, что это Зиммерман написал библиотеку, что во всех режимах имеет правильное округление (кроме комплексных чисел, complex.h), т.е. даже не 0.500. Точнее. Как сказано в статье

      "Note that correct rounding implies an entry 0.5, but the converse is false: when the “round-to-even” rule applies, if the library returns the “odd” value, the error will still be 0.5 ulp." https://gitlab.inria.fr/core-math/core-math/-/commit/6ee58266fe4543511a3444ac4dff38318291cfac
      Часть этого кода Зиммерман апстримит в llvm.


  1. mynameco
    12.02.2025 13:23

    вспоминаю тред, где обсуждалось что простые операции недетерменированные. что один бинарный код, на одной машине, может давать разные результаты. это зависит к какому ядру привязан поток. у каждого потока свой доступ к fpu или общий может быть занят и выполняется на своем. а у общего fpu могут быть более длинные регистры.

    а сейчас со смартфонами, с разношерстными ядрами в одной системе, может быть что угодно. наверняка fpu общий, может отключаться по режимам энергопотребления.

    т.е. это не проблема даже кода. это проблема железа.


    1. mynameco
      12.02.2025 13:23

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

      в одной был явный фиксед поинт класс. на основе инта. во второй игре был просто инт или шорт, без оберток и считалось от мин до макс числа. почти как в марио на денди.


  1. Refridgerator
    12.02.2025 13:23

    Оказалось, что у большинства современных 64-битных процессоров есть достаточно мощные модули для выполнения операций с 64-битными числами с плавающей запятой

    И даже с 80-битной есть (у FPU), использую её в критичных к точности вычислениях.


    1. Panzerschrek
      12.02.2025 13:23

      В своё время этот 80-битный FPU был источником проблем в точности алгоритмов - компилятор мог иногда эти 80 бит до 64 или 32 округлять, а иногда не делал этого. Из-за этого результаты вычислений могли быть разными в зависимости от версии компилятора и его настроек. Также код ломался при переходе на float вычисления на основе SSE инструкций.


      1. Refridgerator
        12.02.2025 13:23

        Я пишу сразу на ассемблере, поэтому проблем таких не имею. Собственно, даже для float/double писать на ассемблере - единственная гарантия детерминированности вычислений.


        1. redfox0
          12.02.2025 13:23

          Я всё-таки поинтеерсуюсь: а почему не используете числа с фиксированной точкой?


          1. Refridgerator
            12.02.2025 13:23

            Ну потому что в числах с фиксированной точкой есть одна проблема - по сути это тот же int. В моих задачах и такое встречалось, и такое я тоже использовал. Но конкретно в этой статье - речь о плавающей точке.


  1. Panzerschrek
    12.02.2025 13:23

    Недетерминированность float вычислений - это не новость. Посему рекомендуется там, где нужна детерминированность, использовать только целочисленные вычисления, ну или эквивалентные им вычисления с фиксированной точкой. Ну кроме случаев, когда целевое оборудование гарантирует детерминированность float операций.


    1. mentin
      12.02.2025 13:23

      Видел доклад на последнем cppcon, кросс-платформенный детерминизм для плавающей точки. Но понятно это значит своя библиотека для всего

      https://github.com/CppCon/CppCon2024/blob/main/Presentations/Cross-Platform_Floating-Point_Determinism_Out_of_the_Box.pdf


  1. Nick0las
    12.02.2025 13:23

    Интересная заметка. Но в глубину главный вопрос не раскрыт. В чем причина разного округления разных реализаций если копнуть до инструкций процессора? По идее и та и другая использует для вычислений FPU, который вычисляет синус одной инструкцией и конвертирует FP32/FP64/FP80 при загрузке в стек тоже единобразно. Или кто-то использует векторные иснтрукции? Или настройки округления в FPU меняются неконтролируемо. Или вообще вместо аппаратного FPU используется эмуляция?

    И это вы еще не сравнивали разные реализации FPU на разных платформах, там, возможно, тоже чудеса будут.


    1. unreal_undead2
      12.02.2025 13:23

      Или кто-то использует векторные иснтрукции? 

      Весь тяжёлый код, который бегает на больших кластерах - однозначно на векторных инструкциях и всяческие синусы реализуются вручную на смеси таблиц и Ньютона-Рафсона - так быстрее, чем на x87.


  1. vadimr
    12.02.2025 13:23

    Тут надо бы начать с того, что само железо реальных процессоров не всегда обеспечивает погрешность в 0.5 ULP. Так что цель благородна, но недостижима.

    Что касается практического примера автора, то, если у него поведение программы зависит от содержимого младшего бита вещественных чисел, то он просто неправильно выбрал разрядность или его алгоритм ошибочен. Хорошо, допустим, он достиг бы идеальной точности в 32-м бите. Но 33-й бит по-прежнему неточен, просто потому, что он не кодируется. Он просто чуть-чуть уменьшил относительную погрешность, а не устранил проблему в принципе.


    1. unC0Rr
      12.02.2025 13:23

      если у него поведение программы зависит от содержимого младшего бита вещественных чисел, то он просто неправильно выбрал разрядность или его алгоритм ошибочен

      Зависимость от младшего бита неизбежна, как только в программе появляется ветвление, основанное на сравнении больше-меньше двух чисел.


      1. unreal_undead2
        12.02.2025 13:23

        Если из за того, что мы каком то проценте случаев из за младшего бита пойдём по другой ветке, существенно (больше заданного порога) меняется выходной результат - надо менять точность промежуточных данных (или сам алгоритм).


        1. iShrimp
          12.02.2025 13:23

          Главное, чтобы результат sin/cos не выходил за пределы [-1; 1].

          В остальном, погрешность в 1-2 бита программиста волновать не должна.


          1. sci_nov
            12.02.2025 13:23

            Вахах)


          1. unreal_undead2
            12.02.2025 13:23

            Если на входе комплексный аргумент - может куда угодно выйти )


        1. unC0Rr
          12.02.2025 13:23

          Если речь об игре, то от результата сравнения может зависеть, есть коллизия между объектами или нет, из чего могут быть очень далеко идущие последствия.


          1. unreal_undead2
            12.02.2025 13:23

            Ну какие могут быть последствия в игре, это ж не управляющая система на спутнике или ядерном реакторе. Максимум падение прибыли разработчиков, да и то не факт, может с ошибкой веселее получится )


        1. grigorym
          12.02.2025 13:23

          Есть вычислительные алгоритмы неустойчивые, потому что неустойчив сам процесс, который они моделируют (детерминированный хаос). У нас, например, это моделирование движения жидкости, где на каждом временном шаге решаются системы уравнений, включающие данные с предыдущего шага. Изменения в младшем бите - мизерные, но от шага к шагу они накапливаются, и результат получается другой.


      1. vadimr
        12.02.2025 13:23

        Однако для вещественных чисел существует понятие погрешности. И если в программе принимается какое-то решение на основании различия чисел, не превышающего их погрешность, то это ошибочная программа. Такие числа следует считать равными.

        Сравнение вещественных чисел должно быть примерно таким:

        if ((a-b) > eps) ...

        вместо

        if (a > b) ...


        1. unreal_undead2
          12.02.2025 13:23

          Первое условие тоже может меняться при изменении младшего бита в a или b.


          1. vadimr
            12.02.2025 13:23

            Формально может, однако eps же выбирается не просто так, а чтобы быть каким-то порогом для реально возможных шумов.


            1. unreal_undead2
              12.02.2025 13:23

              Ни разу не понимаю, чем a > b + eps лучше a > b - просто поменяли условие (и почему не на a > b - eps?). Как вариант - оставить просто a > b, но накапливать статистику случаев abs (a - b ) < eps и если их заметное количество - работать дальше над численной устойчивостью (скажем, увеличивать разрядность a и b, уменьшая eps).


              1. vadimr
                12.02.2025 13:23

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


                1. unreal_undead2
                  12.02.2025 13:23

                  Вопрос, насколько сильно на конкретный алгоритм повлияет выполнение другой ветки при близких a и b. Тут сложно рассуждать абстрактно, каждый случай надо анализировать отдельно.


                  1. SADKO
                    12.02.2025 13:23

                    Всё


                  1. SADKO
                    12.02.2025 13:23

                    Да, есть варианты, но в общем и целом всё гораздо проще ибо FP само по себе является источником погрешностей\шумов, но даже в ситуации когда их вроде-бы не набегает, на результат сравнения они действуют!
                    А сравниваем мы не из любви к искусству, а ради каких-то меркантильных целей. Вот торговый робот, что-то должен сделать если результат одной цепи вычислений, больше результата другой, и НЕ СДЕЛАЕТ. Тк тот будет меньше на какую-то не мыслимую, но критичную для прямого сравнения малую долю.
                    Целочисленные системы отработают быстрее и получат свой гешефт.


                    1. unreal_undead2
                      12.02.2025 13:23

                      Если надо оперировать с величинами, отличающимися на десятки порядков (но при этом при суммировании значений порядка 1e-10 может набежать 1e10) - целые числа и fixed point мало помогут.


    1. Ivan22
      12.02.2025 13:23

       само железо реальных процессоров не всегда обеспечивает погрешность

      мало того, ученые говорят сама реальность вселенной недерминированна и некоторые параметры в квантовом мире в принципе нельзя измерить точно


  1. QuarkDoe
    12.02.2025 13:23

    Для GCC и CLang'а кто-то подобное делал?


  1. ErmIg
    12.02.2025 13:23

    Полная идентичность операций с плавающей запятой не достижима уже на уровне A×B+C != FMA(A,B,C). Т.е. код с AVX2 и без него даёт разную погрешность.