Comments 44
Если нужна точность, то зачем работать с числами типа float?
Для этого есть double и long double.)
Описанные проблемы всплывут и с этими типами, особой разницы нет. Но вообще позиция автора близка к параноидальной
только в самых ограниченных ситуациях при вычислениях суммы мы знаем, что все входные числа являются конечными и что их сумма не может привести к переполнению.
На мой взгляд это скорее common case.
Если вы используете в качестве целевой платформы для своего языка LLVM, то избегайте флагов fast-math
nnanиninf
При расчёте параметров ядерного реактора или ракетного двигателя - однозначно. В игрушках и подобном софте можно и пожертвовать корректной обраткой nan/inf.
с double проблем меньше )
Кстати если точность является критичной, то имеется библиотеки для вычислений с плавающей запятой многократной точности, например, GNU MPFR Library.
Я думаю float это в статье для простоты демонстрации. На практике конечно [long] double, но их тоже не хватает. MPFR будет практически всегда медленнее, чем алгоритмы для конкретного применения, вроде алгоритма Кэхэна для суммы.
Всё зависит от конкретной цели, где это применяется. Где нужно быстродействие, там проигрывает точность, и наоборот. Пока широко не внедрят специализированные АЛУ для работы с числами произвольной точности, говорить об одновременном быстродействии и точности машинной арифметики считаю не имеет смысла.
Вот статья как раз показывает, что такой подход для конкретного (но очень популярного) применения - не правилен, АЛУ не нужно, числа произвольной точности в данных применениях не нужны. Есть алгоритмы, позволяющие рассчитать максимальную погрешность, и выбрать самый быстрый алгоритм обеспечивающий нужную вам точность. И этот алгоритм будет гораздо быстрее, чем код "тупо" использующий MPFR или другие библиотеки чисел с расширенной точностью (которыми тоже надо правильно пользоваться, чтобы правильно рассчитать точность чисел, дающую нужную точность результата).
В принципе соглашусь с вами насчёт выбора алгоритма Кэхэна для суммы. Однако всё же следует ориентироваться с выбором алгоритма в зависимости от предъявляемых требований к задаче, так как возможен вариант, что стандартные типы чисел с плавающей запятой изначально не смогут удовлетворить точность вычисления.
На мой взгляд это скорее common case.
Угу. Если у вас во входном массиве есть хоть один NaN или бесконечность (что технически одно и то же, inf это одно из значений NaN), то, как ни суммируй, на выходе будет NaN.
Если уж прям надо обрабатывать случай с NaN, то можно делать проверку экспоненты каждого числа битовой маской и сразу возвращать это число в качестве суммы, если все биты экспоненты единичны. Поскольку для AND нужно целочисленное АЛУ, на современных процессорах никакого замедления эта проверка не даст.
И всё же с бесконечными (±inf) не всё так однозначно. В качестве слагаемого он даст неопределённость (NaN) только с inf противоположного знака, либо с NaN.
Если же inf окажется в знаменателе, то запросто приведёт результат к обычному ±0.
В статье, как мне кажется, очень неплохо раскрыто, что практически даром можно заложить корректное поведение кода в пограничных случаях
Ещё, конечно, есть варианты с предварительной сортировкой массива. В идеале по возрастанию абсолютного значения, но можно и находить "середину" массива, отсортированного по возрастанию или убыванию, бисекцией и дальше считать в обе стороны.
O(n log n) для сортировки это не так дорого, если например ограничиться блоками в 256 или там 4096 элементов (чтобы влезало в кэш L1). Разумеется, это нужно делать рекурсивно, т.е. на уровне выше сортировать значения сумм по 256 элементов. Точность при этом будет на уровне accurate алгоритмов, а скорость - на уровне блочных.
Работал в финтехе одно время, вот там наверное самые популярные грабли у начинающих разработчиков - это представление сумм денег в виде числа с плавающей запятой, а потом огребание вот этих самых проблем. Потому что складывать надо много, например все транзакции платежной системы за недели/месяцы - это очень много как самих операций сложения, так и потенциально большие числа сами по себе.
Но там это обходят по-своему, используют целочисленную арифметику и ведут расчеты в самых маленьких денежных единицах, например копейках (ну или я в Нигерии работал, там были наиры и кобо). То есть хранят 20 рублей 30 копеек как целое число 2030, а не как вещественное 20.3, и никаких проблем.
В аналитике и моделировании (когда дело не касается реальных сумм на счетах) FP там всё таки уместен - скажем, классическая модель Блэка-Шоулза на float/double считается (хотя я сталкивался только с бенчмарками, а не реальным применением). А с выбором базы в целочисленных расчётах тоже надо действовать аккуратно, особенно если обрабатываются разные валюты, акции и т.д. - навскидку в официальных курсах ЦБ четыре знака после запятой.
А как насчёт вот этого?
Зачем воздвигать самому себе препятствия, используя там где нужна точность числа с плавающей запятой, вместо чисел с фиксированной запятой, а затем героически их преодолевать?
Нууу... числа с фиксированной прихрдитмя придумывать самому
Если обрабатываемые числа отличаются на десятки порядков, никуда не денешься.
Попробуйте прибавить к 1E30 единицу в формате с плавающей запятой.
А числа с фиксированной точкой могут оперировать и с сотнями порядков.
Одна единица на результат не влияет, если их много - см. текст статьи )
А числа с фиксированной точкой могут оперировать и с сотнями порядков.
Что при этом с производительностью?
Одна единица на результат не влияет, если их много - см. текст статьи )
Если речь о 1E30, то даже миллион единиц не повлияет. Так что если нужна точность - необходимо отказаться от плавающей запятой.
Что при этом с производительностью?
Если сравнивать просто SUM() в PostgreSQL, то на 30-35% медленней. Но если вместо SUM() попробовать алгоритмы из текста статьи, то фиксированная запятая сразу выигрывает с отрывом.
Проверять на SIMD, если честно, в лом. Но есть обоснованные подозрения, что выиграет фиксированная запятая. За счет того, что сложение с фиксированной запятой в AVX-512 всегда один такт, а с плавающей - нет, так как требует предварительной нормализации.
Если речь о 1E30, то даже миллион единиц не повлияет.
А вот 1E30 - повлияют, а это не такой уж большой размер данных, на одну ноду влезет.
Если сравнивать просто SUM() в PostgreSQL
Странный референс, я бы скорее на BLAS (или его аналоги и оптимизированные реализации от вендоров железа) смотрел.
За счет того, что сложение с фиксированной запятой в AVX-512 всегда один такт, а с плавающей - нет
Спецификации процессоров с вами не согласны. И не очень понимаю, как за один такт обработать операцию с фиксированной запятой с сотнями порядков.
А вот 1E30 - повлияют, а это не такой уж большой размер данных, на одну ноду влезет.
У Вас с арифметикой проблемы. 4E30 байт - это 4 кветтабайт. При том, что объем всей информации в мире оценивается, примерно, в 100 зеттабайт 1E23. В 10 миллионов раз меньше.
Странный референс, я бы скорее на BLAS (или его аналоги и оптимизированные реализации от вендоров железа) смотрел.
А я Вас не понимаю. BLAS - это просто библиотека. Как Вы в ней будете хранить большие массивы информации? Особенно, если Вы на кветтабайты замахиваетесь )))
И не очень понимаю, как за один такт обработать операцию с фиксированной запятой с сотнями порядков.
512-битное число с фиксированной запятой, вполне умещающееся в регистр ZMM, эквивалентно, примерно, 153 значащим десятичным цифрам. Гугол легко влезет. А больше и не надо, так как это намного больше количества элементарных частиц в видимой части вселенной.
4E30 байт - это 4 кветтабайт
Sorry, по основанию 2 считал )
BLAS - это просто библиотека. Как Вы в ней будете хранить
Мы вроде про стоимость вычислений говорим, а не про хранение.
512-битное число с фиксированной запятой, вполне умещающееся в регистр ZMM
То есть целый регистр (и соответствующее число операций) тратится на одно число вместо восьми даблов - причём явных операций для сложения 512-битных чисел нет, придётся ещё вручную с переносами разбираться.
А больше и не надо, так как это намного больше количества элементарных частиц в видимой части вселенной.
Если у нас исходно целые числа - да, если какой нибудь дифур считаем, разброс порядков никак не связан с числом чего бы то ни было.
Мы вроде про стоимость вычислений говорим, а не про хранение.
Одно без другого не бывает. Особенно с учетом того, что в статье рассматриваются методы не последовательного суммирования.
Нравится Вам это или нет, но числа откуда-то берутся. И в большинстве случаев - именно из базы данных. Да и мне проще было сгенерировать таблички из случайных чисел и посчитать время по ним.
Если у нас исходно целые числа - да, если какой нибудь дифур считаем, разброс порядков никак не связан с числом чего бы то ни было.
В любом случае, число больше гугола никакого смысла не имеет.
Нравится Вам это или нет, но числа откуда-то берутся. И в большинстве случаев - именно из базы данных
Это могут промежуточные результаты вычислений в памяти (тот же дифур на большой сетке).
Я потому и написал, что "в большинстве случаев", а не всегда.
Ну привели Вы исключение, когда данные могут уместиться в оперативной памяти для последующего применения изложенных в статье методов. Что это меняет?
На таких объемах эти миллисекунды все рано роли не сыграют.
Тут вопрос как на мир смотреть, для меня большие объёмы действительных чисел - это в первую очередь HPC и сложные расчёты.
На таких объемах эти миллисекунды все рано роли не сыграют.
Смотря насколько часто эти вычисления проделываются и сколько раз, миллисекунды могут и в часы сложиться (хотя и просто пробежаться по нескольким сотням гигабайт - уже не миллисекунды).
для меня большие объёмы действительных чисел - это в первую очередь HPC и сложные расчёты
Так по любому миллиарды чисел в оперативку не влезут. А миллионы суммировать - десятки миллисекунд.
Смотря насколько часто эти вычисления проделываются и сколько раз, миллисекунды могут и в часы сложиться
Вы лукавите. Тут опять, либо это параллельное суммирование на разных ядрах, что еще больше мешает держать все данные в оперативке, либо несколько последовательных, где из миллисекунд на суммировании даже 1% вряд ли сложится.
Или Вы пытаетесь подменить тему, переходя от суммирования очень длинного вектора к совсем другим расчетам?
хотя и просто пробежаться по нескольким сотням гигабайт - уже не миллисекунды
Сотни гигабайт - это уже явно вектор не в памяти, а в БД. А вот миллисекунды или нет - зависит от мощности кластера.
А вообще Вы мне напоминаете экскурсовода, утверждающего, что этой окаменелости 200 000 005 лет, потому что пять лет назад ей было 200 миллионов лет. Точные вычисления - это не про дифуры, где исходные данные хорошо, если шесть точных знаков имеют.
Так по любому миллиарды чисел в оперативку не влезут.
Берём, скажем, такую машинку - https://www.alcf.anl.gov/aurora - 512 Gb обычной DRAM на узле плюс ещё сколько то HBM, вполне хватит для миллиардов чисел.
переходя от суммирования очень длинного вектора к совсем другим расчетам?
Понятно, что суммирование вектора - только одна из промежуточных подзадач в сложном расчёте (скажем, получить общую энергию на сетке после нескольких шагов для проверки, что считаем правильно).
Точные вычисления - это не про дифуры, где исходные данные хорошо, если шесть точных знаков имеют.
Так речь не об абсолютной точности, а чтобы не получить ошибку в разы, которую может дать сложение в лоб.
Берём, скажем, такую машинку
И вычитаем из Вашей заработной платы её амортизацию, так как не хватило профессионализма воспользоваться СУБД )))
одна из промежуточных подзадач в сложном расчёте
Я и говорю, подмена темы.
Еще раз, если окаменелости пять лет назад было 200 млн. лет, то сейчас ей 200 000 005 лет? Тоже ведь единицу каждый год добавляем? )))
не получить ошибку в разы
Просто складывайте миллиард единиц (9 знаков) в f32, точность которого ограничена 7 десятичными знаками.
так как не хватило профессионализма воспользоваться СУБД
Для расчёта, скажем, термоядерной реакции?
Тоже ведь единицу каждый год добавляем?
Если добавили пять - то так и будет 2e8. Если моделируем динамику на миллиард лет - таки надо будет сложить единички )
Просто складывайте миллиард единиц (9 знаков) в f32, точность которого ограничена 7 десятичными знаками.
Подозреваю, что имелось в виду "не cкладывайте" - но и с даблом накопленная погрешность может быть больше, чем устраивает.
Для расчёта, скажем, термоядерной реакции?
Есть другие варианты? CERN именно в БД результаты и загоняет. Потому что обрабатывать их в реальном времени невозможно физически.
погрешность может быть больше
Для того и есть числа с фиксированной запятой, точность которых ограничена лишь объемом памяти.
CERN именно в БД результаты и загоняет. Потому что обрабатывать их в реальном времени невозможно физически.
Это обработка результатов измерений, я про теоретические расчёты.
Для того и есть числа с фиксированной запятой, точность которых ограничена лишь объемом памяти.
При этом и сложность той же операции сложения двух чисел будет расти пропорционально. Возможно, в каких то ситуациях тратить гигабайты на одно число имеет смысл, но чаще есть компромисс между точностью и скоростью, так что в зависимости от ситуации может просто хватить float или double, можно запользовать float128, можно явно отслеживать ошибки округления в алгоритме (как при сложении Кэхена). Фиксированная точка тоже имеет право на жизнь, но про широкое использование в HPC не слышал.
Это обработка результатов измерений, я про теоретические расчёты.
которые или возникают на основании результатов измерений или проверяются ими.
При этом и сложность той же операции сложения двух чисел будет расти пропорционально.
Линейно, в зависимости от востребованной точности. А вот сложность алгоритмов сложения в статье вовсе даже не линейна. Плюс еще трудозатраты, которые почти нулевые при вызове SUM() и могут быть весьма велики в других алгоритмах.
Я же провел эксперимент. Повторите. Создайте таблички с миллиардом значений float и decimal в PostgreSQL со случайными величинами. А потом попробуйте написать запрос, которых суммирует все эти значения быстрее, чем встроенный SUM() по decimal и не накапливает погрешность. Вот тогда можно будет продолжить дискуссию.
про широкое использование в HPC не слышал
Потому что 15 десятичных знаков достаточно в большинстве случаев, а 34/71 десятичных знаков в f128/f256 достаточно в подавляющем большинстве случаев.
Другой способ заключается в хранении большого буфера целых чисел...
Этот способ требует времени 𝑂(𝑛), но использует большой, хотя и константный объём памяти (≈≈ 1 КБ в случае f32, ≈≈ 16 КБ в случае f64).
Это неправильная оценка. В f32 8 бит экспонента (порядок) и 23 мантисса. Итого достаточно 279 бит. Не то чтобы много, но никаких килобайт не надо. Это 9 32-битных чисел, 5 64-битных или 3 128-битных xmm регистра. В f64 11 бит порядок и 52 мантисса. 66 32-битных, 33 64-битных, 17 128-битных чисел.Опять же - ужас, но не ужас-ужас.
Другой вопрос, что по ним распихивать слагаемые неудобно (навскидку - 5-10 тактов запросто уйдёт).
Откуда вы взяли 279 бит? Экспонента float32 может иметь 256 значений, по 4 байта на каждое значение это ровно 1 КиБ (ну ок, значения 0 и 255 вырожденные, но это особо не влияет)
Смотрите. Рассмотрим простой-простой "uFP8" тип беззнаковый и 8-битный. 4 бита на порядок, 4 на экспоненту. Забудем о NaN, Inf и нормализацию и уникальность нуля. Будем считать, что экспонента всегда положительная. Будем считать, что старший бит мантиссы соответствует 1*2^exp, где exp это порядок. Упрощений много, но на самом деле они почти не влияют на конечную картину. Т.е. у нас есть 0 (0000-0000, где первые 4 бита порядок, потом мантисса) и числа от 1/8 (0000-0001) до 15*2^15=491520 (1111-1111).
При этом из-за отсутствия нормализации некоторые числа имеют несколько представлений. Например, 4 может быть представлено как 0010-1000 или как 0011-0100 или как 0100-0010 или как 0101-0001.
Теперь будем раскладывать наше uFP8 число, например 1280 (1010-1010) на 4-битные непересекающиеся компоненты:
Для примера, самое маленькое число:
0000-0001 xxxxxxxxxxxxxxx0001 - это 1/8, на местах x неявные нулевые биты
Самое большое число:
1111-1111 1111xxxxxxxxxxxxxxx - это 15*2^15=491520, на местах x неявные нулевые биты
Наше число:
1010-1010 xxxxx1010xxxxxxxxxx
порядок мантисса
0000 0000
0100 0000
1000 1000
1100 0010
10000 x000 - тут на месте x всегда неявный 0
Таким образом мы каждое наше исходное число можем разложить в 4 числа по 4 бита и одно по 3, т.е. в 19 бит.
Технически (для управления переполнениями при сложении), наверное чуть удобнее раскладывать в 4 пятибитовых числа и 1 четырёхбитовое - 24 бита суммарно. Но это непринципиально. Ну и да, раскладка сожрёт несколько тактов.
Аналогично, но с учетом других размеров, знака, NaN, Inf и нормализации можно раскладывать и FP32/FP64.
Нельзя ли просто привести к целочисленному типу и уже работать с ним?
Статья полезна для ознакомления, но на практике применить особо негде. Заказчик не поймёт таких трат. Поэтому для денег используем decimal (который из шарпа), его аналоги или целые числа, для прочих - float
Укрощаем суммы с плавающей запятой