Comments 67
Фундаментально!
похоже, мы можем обеспечить соответствие стандарту IEEE 754 в 32-битных реализациях на современных машинах без существенного ущерба производительности.
Похоже, тестировался невекторизованный код.
Статья 5ти летней давности. Вот что-то не такое старое https://members.loria.fr/PZimmermann/papers/accuracy.pdf
Вы главное забыли упомянуть, в этой статье упомянуто, что это Зиммерман написал библиотеку, что во всех режимах имеет правильное округление (кроме комплексных чисел, 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.
вспоминаю тред, где обсуждалось что простые операции недетерменированные. что один бинарный код, на одной машине, может давать разные результаты. это зависит к какому ядру привязан поток. у каждого потока свой доступ к fpu или общий может быть занят и выполняется на своем. а у общего fpu могут быть более длинные регистры.
а сейчас со смартфонами, с разношерстными ядрами в одной системе, может быть что угодно. наверняка fpu общий, может отключаться по режимам энергопотребления.
т.е. это не проблема даже кода. это проблема железа.
реверс инженерил две популярные мобильные игры давно. асинхронный мультиплеер, там нужны детерминированные симуляции. чтобы сервак мог проверить, мог ли ты так, разрушить замок и побить врага.
в одной был явный фиксед поинт класс. на основе инта. во второй игре был просто инт или шорт, без оберток и считалось от мин до макс числа. почти как в марио на денди.
Вот да. Играл когда-то в 3D шутер, уперся в тупик - нужно было перепрыгнуть далеко, к тому же с набором высоты - никак не получалось. На компьютере товарища получилось с первого раза. Проверил на своем - снова никак. Тогда сделал вывод, что в условиях дефицита производительности игра пропускала итерации расчета координат, чтобы не было критического провала fps. А может все же не пропуски, а упрощенный алгоритм расчета.
что один бинарный код, на одной машине, может давать разные результаты. это зависит к какому ядру привязан поток.
Вот в это крайне слабо верится. Можете найти ссылку?
https://members.loria.fr/PZimmermann/papers/accuracy.pdf
However, it was observed that the Intel Math Library gives better results on AMD hardware than on Intel hardware for acosh, asin, asinh and atan2; a possible explanation is that these functions use the rsqrt instruction, which is known to be more accurate on AMD hardware https://github.com/jeff-arnold/math_routines.git
Кроме того: https://members.loria.fr/PZimmermann/talks/core-math-arith2022.pdf
страница 4
И где там про зависимость от привязки потока к ядру на конкретной модели процессора?
Ну да, на разных процессорах операции типа не-совсем-элементарных функций (FSIN) или заведомо приблизительные (RSQRT*) могут выполняться по-разному. Это известно.
Но чтобы это было на разных ядрах одного процессора - я ещё такого не слышал и ваша ссылка не даёт. Даже с современным разделением на P-cores и E-cores они должны работать в таких вопросах абсолютно одинаково. Так что не засчитывается.
Оказалось, что у большинства современных 64-битных процессоров есть достаточно мощные модули для выполнения операций с 64-битными числами с плавающей запятой
И даже с 80-битной есть (у FPU), использую её в критичных к точности вычислениях.
В своё время этот 80-битный FPU был источником проблем в точности алгоритмов - компилятор мог иногда эти 80 бит до 64 или 32 округлять, а иногда не делал этого. Из-за этого результаты вычислений могли быть разными в зависимости от версии компилятора и его настроек. Также код ломался при переходе на float вычисления на основе SSE инструкций.
Я пишу сразу на ассемблере, поэтому проблем таких не имею. Собственно, даже для float/double писать на ассемблере - единственная гарантия детерминированности вычислений.
Я всё-таки поинтеерсуюсь: а почему не используете числа с фиксированной точкой?
В случае морфингового процессора (Crusoe, который мог менять пол архитектуру) тоже?
К примеру, прилетело обновление винды, обновился микрокод процессора и появились некоторые особенности обработки машинных кодов.
Или злобные хакеры прошлись хаммером по ОЗУ.
*фантазии, навеянные фразой "гарантия детерминированности вычислений"
"компилятор мог иногда эти 80 бит до 64 или 32 округлять, а иногда не делал этого"
Вы наверное все же имеете ввиду bug 323
Недетерминированность float вычислений - это не новость. Посему рекомендуется там, где нужна детерминированность, использовать только целочисленные вычисления, ну или эквивалентные им вычисления с фиксированной точкой. Ну кроме случаев, когда целевое оборудование гарантирует детерминированность float операций.
Интересная заметка. Но в глубину главный вопрос не раскрыт. В чем причина разного округления разных реализаций если копнуть до инструкций процессора? По идее и та и другая использует для вычислений FPU, который вычисляет синус одной инструкцией и конвертирует FP32/FP64/FP80 при загрузке в стек тоже единобразно. Или кто-то использует векторные иснтрукции? Или настройки округления в FPU меняются неконтролируемо. Или вообще вместо аппаратного FPU используется эмуляция?
И это вы еще не сравнивали разные реализации FPU на разных платформах, там, возможно, тоже чудеса будут.
Или кто-то использует векторные иснтрукции?
Весь тяжёлый код, который бегает на больших кластерах - однозначно на векторных инструкциях и всяческие синусы реализуются вручную на смеси таблиц и Ньютона-Рафсона - так быстрее, чем на x87.
По идее и та и другая использует для вычислений FPU, который вычисляет синус одной инструкцией и конвертирует FP32/FP64/FP80 при загрузке в стек тоже единобразно.
Вы таки сильно отстали. Microsoft в 64-битном режиме FPU вообще не допускает, только SSE. В Linux, FreeBSD и т.д. возможно, но с использованием long double или явного ключа (GCC) -mfpmath=387, который по умолчанию в 64-битном режиме не включается. В 32 битах, да, было FPU, если не сказано перейти на SSE (-march=i686 как минимум плюс -mfpmath=sse). При использовании SSE функции типа синуса реализуются библиотекой, а не процессором.
Называть SSE "векторными инструкциями" некорректно, там специально выделена группа скалярных на тех же регистрах. Но если отвлечься от этого - да, "кто-то", то есть сейчас чуть менее чем все, использует именно их.
Вопросы других платформ не рассматриваю, там своя специфика у каждой.
И это вы еще не сравнивали разные реализации FPU на разных платформах, там, возможно, тоже чудеса будут.
Конечно, будут. Реализации за пределами 5 операций арифметики, точного квадратного корня и конверсий типа отличаются уже у Intel и AMD (по ссылке как раз то, что активно используется в играх). А конверсии - и между ISA. Например, если число во float не влезает в int, то приведение к целому даёт ближайшее крайнее значение у ARM, но специальный маркер на x86. Таких примеров много. Но статья не об этом.
По идее и та и другая использует для вычислений FPU, который вычисляет синус одной инструкцией и конвертирует FP32/FP64/FP80 при загрузке в стек тоже единобразно.
Вы таки сильно отстали. Microsoft в 64-битном режиме FPU вообще не допускает, только SSE. В Linux, FreeBSD и т.д. возможно, но с использованием long double или явного ключа (GCC) -mfpmath=387, который по умолчанию в 64-битном режиме не включается. В 32 битах, да, было FPU, если не сказано перейти на SSE (-march=i686 как минимум плюс -mfpmath=sse). При использовании SSE функции типа синуса реализуются библиотекой, а не процессором.
Называть SSE "векторными инструкциями" некорректно, там специально выделена группа скалярных на тех же регистрах. Но если отвлечься от этого - да, "кто-то", то есть сейчас чуть менее чем все, использует именно их.
Вопросы других платформ не рассматриваю, там своя специфика у каждой.
И это вы еще не сравнивали разные реализации FPU на разных платформах, там, возможно, тоже чудеса будут.
Конечно, будут. Реализации за пределами 5 операций арифметики, точного квадратного корня и конверсий типа отличаются уже у Intel и AMD (по ссылке как раз то, что активно используется в играх). А конверсии - и между ISA. Например, если число во float не влезает в int, то приведение к целому даёт ближайшее крайнее значение у ARM, но специальный маркер на x86. Таких примеров много. Но статья не об этом.
В чем причина разного округления разных реализаций если копнуть до инструкций процессора?
Во-первых, разные процессоры, причём разных архитектур. Во-вторых, идеальная точность это очень сложно в разработке (требуется серьёзная наука) и дорого в реализации. Рядом вот ссылаются на новые работы в этой области - но до популярных библиотек результаты ещё не дошли.
Называть SSE "векторными инструкциями" некорректно, там специально выделена группа скалярных на тех же регистрах.
Да, если ничего специально не предпринимать, будут работать скалярные инструкции с одним элементом вектора. Но если код вычислительно сложный и реально нужный, скорее всего над его оптимизацией работали знающие люди (возможно, инженеры того же Intel) - и тогда там реальные вектора, и не SSE, а AVX или AVX512.
и тогда там реальные вектора, и не SSE, а AVX или AVX512.
Я понимаю вашу мысль и частично согласен, но:
1) Если код не содержит массивно-параллельных операций, то векторизовать нечего. Или может быть, что векторизация невозможна из-за сложности кода (в основном, чрезмерном количестве развилок в логике).
2) Если код содержит много массивно-параллельных операций, то его выгоднее выполнять на GPGPU. Лучше выгнуться и потратиться на оптимизацию под такое применение, но пойти туда, где будет не AVX512, а условный эквивалент AVX65536:)) и инженеры будут не от Intel, а от NVidia.
Для SSE/AVX/etc. в процессоре остаются ниши или промежуточные между этими, или там, где торможение от загрузки данных в GPGPU и обратно дороже, чем исполнение на основном процессоре.
Проблема векторности в процессоре на сейчас именно в этом. И вот для "элементарных функций" из исходной статьи переложить на векторизуемый вид - это отдельная научная работа.
И вот для "элементарных функций" из исходной статьи переложить на векторизуемый вид - это отдельная научная работа.
Да обычная инженерная работа, в SVML всё давно есть, интеловский компилятор может её автоматом запользовать.
Тут надо бы начать с того, что само железо реальных процессоров не всегда обеспечивает погрешность в 0.5 ULP. Так что цель благородна, но недостижима.
Что касается практического примера автора, то, если у него поведение программы зависит от содержимого младшего бита вещественных чисел, то он просто неправильно выбрал разрядность или его алгоритм ошибочен. Хорошо, допустим, он достиг бы идеальной точности в 32-м бите. Но 33-й бит по-прежнему неточен, просто потому, что он не кодируется. Он просто чуть-чуть уменьшил относительную погрешность, а не устранил проблему в принципе.
если у него поведение программы зависит от содержимого младшего бита вещественных чисел, то он просто неправильно выбрал разрядность или его алгоритм ошибочен
Зависимость от младшего бита неизбежна, как только в программе появляется ветвление, основанное на сравнении больше-меньше двух чисел.
Если из за того, что мы каком то проценте случаев из за младшего бита пойдём по другой ветке, существенно (больше заданного порога) меняется выходной результат - надо менять точность промежуточных данных (или сам алгоритм).
Главное, чтобы результат sin/cos не выходил за пределы [-1; 1].
В остальном, погрешность в 1-2 бита программиста волновать не должна.
Если речь об игре, то от результата сравнения может зависеть, есть коллизия между объектами или нет, из чего могут быть очень далеко идущие последствия.
Есть вычислительные алгоритмы неустойчивые, потому что неустойчив сам процесс, который они моделируют (детерминированный хаос). У нас, например, это моделирование движения жидкости, где на каждом временном шаге решаются системы уравнений, включающие данные с предыдущего шага. Изменения в младшем бите - мизерные, но от шага к шагу они накапливаются, и результат получается другой.
В таком случае стабильный результат на микроуровне (скорость/давление в конкретном узле сетки) в принципе нельзя стабильно получить, максимум - усреднённые занчение по макрообластям, и они от младшего бита зависеть не должны.
Вопрос, что значит "стабильно"? Устойчиво к малым изменениям входных данных? Да, не получится, потому что это детерминированный хаос, по определению. Но означает ли это, что мы должны тогда интересоваться только усреднёнными значениями? Вовсе нет! Дисперсия детального рисунка очень даже важна. Это как лёд на стекле: да, средняя масса его зависит от температуры и влажности, но рисунок снежинок всё время разный, и он тоже важен. Даже больше скажу: усреднять результаты до средних значений по макрообластям - ошибка.
Однако для вещественных чисел существует понятие погрешности. И если в программе принимается какое-то решение на основании различия чисел, не превышающего их погрешность, то это ошибочная программа. Такие числа следует считать равными.
Сравнение вещественных чисел должно быть примерно таким:
if ((a-b) > eps) ...
вместо
if (a > b) ...
Первое условие тоже может меняться при изменении младшего бита в a или b.
Формально может, однако eps же выбирается не просто так, а чтобы быть каким-то порогом для реально возможных шумов.
Ни разу не понимаю, чем a > b + eps лучше a > b - просто поменяли условие (и почему не на a > b - eps?). Как вариант - оставить просто a > b, но накапливать статистику случаев abs (a - b ) < eps и если их заметное количество - работать дальше над численной устойчивостью (скажем, увеличивать разрядность a и b, уменьшая eps).
В численно устойчивом алгоритме по определению всегда есть такой порог, за который его неустойчивость не выходит. То есть представление меньшего числа не может оказаться слишком уж сильно больше, чеме представление большего числа.
Вопрос, насколько сильно на конкретный алгоритм повлияет выполнение другой ветки при близких a и b. Тут сложно рассуждать абстрактно, каждый случай надо анализировать отдельно.
Всё
Да, есть варианты, но в общем и целом всё гораздо проще ибо FP само по себе является источником погрешностей\шумов, но даже в ситуации когда их вроде-бы не набегает, на результат сравнения они действуют!
А сравниваем мы не из любви к искусству, а ради каких-то меркантильных целей. Вот торговый робот, что-то должен сделать если результат одной цепи вычислений, больше результата другой, и НЕ СДЕЛАЕТ. Тк тот будет меньше на какую-то не мыслимую, но критичную для прямого сравнения малую долю.
Целочисленные системы отработают быстрее и получат свой гешефт.
Может оказаться, что из-за округления b-eps == b (технически). Поэтому, если предполагаем, что числа близкие, берём их разность и уже её изучаем. Но лучше перейти на относительную погрешность. В Питоне, например, math.isclose() делает это в удобном виде.
но накапливать статистику случаев
Так и делают. Но случай разный бывает. Например, итеративный алгоритм даёт значение, близкое к желаемому, но не может идеально приблизиться.
скорее abs(a-b) > eps
Это если надо понять разные числа или близкие. Но логика может быть завязана на то, какое из чисел больше, т .е. просто a > b.
Но логика может быть завязана на то, какое из чисел больше, т .е. просто a > b.
Да, я стормозил. Спасибо, удалил это. Тогда ещё чуть хитрее - у нас есть область точно "да", точно "нет" и область неопределённости (когда 0<=a-b<=eps или abs(a-b)<=eps). В зависимости от задачи надо выбирать для каждого конкретного случая, что мы делаем в случае неопределённости, ну и саму эту область уточнять.
В математическом алгоритме может быть только два варианта, так что случай неопределённости можно разве что залогировать чтобы понять масштаб бедствия.
В математическом алгоритме может быть только два варианта
Ну так то в математическом. Там всё идеально, там и catastrophic cancellation не водится, и какое-нибудь "рекуррентное соотношение Мюллера" (хорошая проба на умение хитрых реализаций) сходится к нужным числам без проблем, и ещё 100500 проблем реального мира не присутствуют.
А вот что делать в промежуточном варианте - масса случаев. Например, мы решаем некое уравнение численно и итеративно, и считаем невязку решения, которая в идеале должна быть ноль, но "in practice they differ" (ц). Шансов, что невязка уйдёт в ноль - в общем случае как раз ноль. А вот дальше вопрос, как описывать критерий остановки итерирования. Например, может быть, что при относительной погрешности до 1e-9 мы делаем ещё максимум 2 итерации, а 1e-6 - 200. Но дальше таки стоп и выдавать уже то, что получилось, с соответствующим плачем.
само железо реальных процессоров не всегда обеспечивает погрешность
мало того, ученые говорят сама реальность вселенной недерминированна и некоторые параметры в квантовом мире в принципе нельзя измерить точно
Для GCC и CLang'а кто-то подобное делал?
Полная идентичность операций с плавающей запятой не достижима уже на уровне A×B+C != FMA(A,B,C). Т.е. код с AVX2 и без него даёт разную погрешность.
Чисто формально, поскольку это всё "_recommended_ operations", то нарушения стандарта нет. Полные реализации могут быть чудовищно дороги и вообще бессмысленны. Например, пусть у нас число 1.0e308, какой вообще смысл брать от него синус? (Вот тут обсуждалось, например.) Но если ты утверждаешь, что соответствуешь стандарту именно в этой части (а не только в арифметической базе), то будешь вынужден считать с почти 1100 битами точности.
FPU x87 декларирован с погрешностью до 1.5 ULP:
With the Pentium processor and later IA-32 processors, the worst case error on transcendental functions is less than 1 ulp when rounding to the nearest (even) and less than 1.5 ulps when rounding in other modes. The functions are guaranteed to be monotonic, with respect to the input operands, throughout the domain supported by the instruction.
Сделать точность выше - растут вычислительные затраты и становится очень серьёзной научной проблемой. Если сделают дёшево - отлично. Если нет - считаю, критической проблемы всё равно нет.
Intel нагло солгала в этой части документа, причём ещё в 2014 исправили они эту ошибку в документации, вы вообще pdf обновляете хоть иногда? Нет, ошибка там в квинтиллионы в fsin... вау. Это вызвало целую пачку скандалов, а вы и не в курсе! Glibc после этого выпилила fsin x86_64 инструкцию. Fsin причём это
микрокодная инструкция, то есть она медленная. https://web.archive.org/web/20141009212414/https://software.intel.com/blogs/2014/10/09/fsin-documentation-improvements-in-the-intel-64-and-ia-32-architectures-software
Errata:
"the worst case error on transcendental functions, when applied to the reduced argument, is less"
https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/
"Сделать точность выше - растут вычислительные затраты и становится очень серьёзной научной проблемой"
Которая уже решена. Glibc недавно получил код от Зиммермана, который быстрее и имеет максимальную точность (то есть не только 0.500 ULP, но в режимах округления выбирает верный из двух 0.500 ULP).
https://members.loria.fr/PZimmermann/talks/core-math-fpbench.pdf
15 патчей для glibc уже в glibc! https://sourceware.org/pipermail/libc-alpha/2025-February/164658.html
AMD скопировали FSIN код полностью.
Подробнее https://notabs.org/fpuaccuracy/
Intel нагло солгала в этой части документа, причём ещё в 2014 исправили они эту ошибку в документации, вы вообще pdf обновляете хоть иногда?
Эта цитата была из документа 325462-084US (combined SDM) версии, обозначенной "June 2024", страница 226. Прежде чем выдвигать обвинение, тем более в агрессивном тоне, проверяйте факты.
(В доке уточнено про sin и прочие, и я особенно подчёркиваю, что я не цитировал это уточнение именно потому, что не вижу смысла в решении подобного запроса.)
Нет, ошибка там в квинтиллионы в fsin... вау. Это вызвало целую пачку скандалов, а вы и не в курсе!
Я в курсе про этот случай и не считаю его сколь-нибудь существенным. Если у вас nextafter(x)-x > π, то я не считаю корректным в принципе требовать конкретный ответ, я бы отдавал NaN. И перестаньте орать, ваши аргументы от этого не станут убедительнее. "Скандалы" по поводу этого разумные люди не выдвигают, это не FDIV bug.
Если же исключить подобные диверсионные случаи и войти в разумные пределы, то и 1.5ULP для таких функций по состоянию на 30 лет назад, и чуть более точное сейчас - это уже неплохо, но, главное, не должно быть критично для правильно написанного алгоритма, который их использует.
5 патчей для glibc уже в glibc! https://sourceware.org/pipermail/libc-alpha/2025-February/164658.html
То есть ещё месяца не прошло, в стабильных версиях ещё нет, реального большого тестирования на широкой массе пользователей не прошло(? - одна частная библиотека, пусть даже лучшая в мире, это слишком мало). В коммитах идёт разное включая "The logic was copied wrong from CORE-MATH.". Конечно, хорошо, что прогресс идёт (да, про эти работы Зиммермана/Циммермана я успел раньше прочитать), но я бы с годик-другой подождал до утверждения, что всё отлично.
Ссылки в вашем комментарии полезные, пусть будут в этой дискуссии, кто-нибудь ещё увидит. Но общий накал лучше бы снизить минимум на пол-порядка.
"(В доке уточнено про sin и прочие, и я особенно подчёркиваю, что я не цитировал это уточнение именно потому, что не вижу смысла в решении подобного запроса.)"
Процитирую за вас:
Сразу в след. параграфе (на самом деле 68 (nearest-even) или 69 (rounded toward zero (truncated)), про это сказано в 8.3.8):
However, for FSIN, FCOS, FSINCOS, and FPTAN which approximate periodic trigonometric functions, the previous
statement about maximum ulp errors is true only when these instructions are applied to reduced argument (see
Section 8.3.8, “Approximation of Pi”). This is due to the fact that only 66 significant bits are retained in the finite
approximation Pi of the number π (3.1415926…), used internally for calculating the reduced argument in FSIN,
FCOS, FSINCOS, and FPTAN. This approximation of π is not always sufficiently accurate for good argument reduction.
"я бы с годик-другой подождал до утверждения, что всё отлично."
Там можно легко перебор сделать single. Ошибок как бы нет.
У пользователей на разных машинах генерировались разные паттерны, то есть это ломало многопользовательский режим нашей игры. Выслеживая ошибку, мы добрались до GPU, имеющих собственные правила, и до погрешности в 1ULP на стороне CPU.
Может проблема здесь больше в том, что управляющие регистры FPU позволяют снижать точность расчётов, а ряд API той же Windows иногда сбивают точность расчётов GPU на самую низкую?
В статье про это ни слова.
Ни одна реализация элементарных функций не соответствует стандарту IEEE 754