Ни одна реализация элементарных функций не соответствует стандарту IEEE 754

Введённый в 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 💻