Однажды передо мной встала задача векторизовать функцию вычисления экспоненты. Неожиданно оказалось, что готового решения не существует. Функции быстрого вычисления экспоненты, использующие векторный код, имеются практически для всех платформ в составе быстрых математических библиотек. Но они, как правило, читают данные из массивов в памяти и возвращают результат обратно в память. А вот такого, чтобы взять данные из регистра и ответ поместить обратно в регистр, не нашлось. Intel, правда, реализовал функцию векторного вычисления экспоненты в своей библиотеке SVML. Microsoft лицензировала эту библиотеку для использования в составе Visual Studio. В этом случае проблем нет. Но если захочется портировать код под GCC, окажется, что SVML в составе стандартных библиотек отсутствует. Пришлось писать свою функцию.
Для начала вспомним как можно вычислить экспоненту. Курс математического анализа подсказывает, что это можно сделать разложением в ряд.
Формула восходит по крайней мере к Эйлеру. И замечательна она тем, что ряд сходится для любого значения x, какой бы большой x не был. Правда, несколько смущает значок бесконечности как верхний предел суммы. Хотелось бы все-таки ограничить число элементов ряда, которые нужно сложить. Для этого нужно ограничить возможные значения x.
Попробуем сначала инженерный подход. В стандарте IEEE 754 число 0 в 32-битном формате представлено следующим образом.
Вполне логично. Но как же быть с правилом, что старший скрытый бит мантиссы всегда равен единице? Нет правил без исключений. Если все биты экспоненты установлены в 0, число считается денормализованным, и скрытый бит полагается равным нулю.
Прибавим единичку к самому младшему биту мантиссы и получим самое маленькое положительное число, представимое в этом формате.
Десятичное представление этого числа 1*10-45, натуральный логарифм его приблизительно равен -103.6.
На другом краю числового кольца лежит значение бесконечности, по стандарту представленное экспонентой, состоящей из единиц, и мантиссой, состоящей из одних нулей.
Вычтем единичку из младшего разряда, и получим самое большое число, на волосок от бесконечности.
Приблизительное десятичное представление этого числа 3.4028235*1038, натуральный логарифм его 88.7.
Интересно, а что будет, если не вычесть, а прибавить единичку к младшему биту бесконечности?
Получится нечисло, Not a Number, NaN. Согласно стандарту IEEE 754, если все биты экспоненты выставлены в 1, а мантисса не равна нулю, то это нечисло. У нечисел есть еще и внутренняя классификация, но эта тема лежит за пределами нашего повествования. Отметим только одно любопытное свойство нечисла: оно не равно ничему. Поэтому самый быстрый способ проверки на нечисло – сравнить значение с самим собой. Если получилось false, мы столкнулись с нечислом.
Наши упражнения с границами представления чисел наводят на простую идею: ограничить x интервалом от -103.0 до +88.0. Идея разумная, поскольку аргумент, выходящий за этот интервал, приведет к выходу результата за рамки представимых 32-битных чисел с плавающей точкой, и поэтому вычисление бессмысленно. Можно сразу сказать, что результат будет либо ноль, если аргумент отрицательный, либо бесконечность, если аргумент положительный. Но, по моей оценке, если мы ограничиваем x столь большим интервалом, нам понадобится больше сотни членов ряда для получения приемлемо точных значений экспоненты, когда факториал заведомо превзойдёт степень аргумента, и ряд начнёт сходиться.
Воспользуемся математической хитростью и представим x в виде целого числа логарифмов двойки плюс некий остаток.
Тогда экспоненту можно вычислить как 2 в степени это самое целое число умноженное на е в степени остаток.
Умножать числа с плавающей точкой на степень двойки мы умеем хорошо. Для этого нужно просто прибавить к экспоненте показатель степени двойки. А g при этом оказывается ограниченным логарифмом двойки, то есть примерно 0,7. Это уже очень приятный интервал аргументов для вычислений, и можно обойтись всего восемью членами ряда.
Представим эту голову ряда виде многочлена.
Вычислим его по школьным формулам.
Помимо векторного параллелизма, школьный подход позволяет также распараллелить вычисления, например, в строках 24 и 25, в строках 26 и 27 и так далее за счет внеочередного выполнения инструкций.
Выполнение программы немного ускоряется за счет использования инструкций комбинированного умножения-сложения (КУС, они же FMA – Fused Multiply Add). При КУС результат умножения прибавляется к сумме нарастающим итогом минуя отсечение дополнительных битов мантиссы в арифметико-логическом устройстве (АЛУ), тем самым уменьшая погрешность округления и экономя время, затрачиваемое на отсечение битов и последующую повторную загрузку в АЛУ.
Еще отметим, что операция умножения на степень двойки настолько важна и распространена, что для нее выделена специальная машинная команда – scalef в строке 36.
Протестируем наш код.
Тест проверяет, чтобы на всем интервале допустимых значений относительная погрешность не превышала заданной величины. Погрешность мы считаем не относительно «математически верного» значения, а относительно результата скалярной библиотечной функции. Ведь при векторизации приходится доказывать, что результат выполнения векторного кода максимально приближен к скалярной версии, уже заслуживающей доверия.
Погрешность более чем в два раза превышает норматив, установленный из теоретического анализа погрешности. Даже КУС не помогло.
Причина накопления погрешности кроется в том, что с ростом степени прибавляемые одночлены становятся все меньше по модулю, и ошибка округления все заметнее.
Алгоритм Кагана, описанный здесь, позволяет учесть и уменьшить возникающую погрешность за счет лишних арифметических операций. Нет ли более дешевых альтернатив?
Компенсировать ошибку округления можно представив многочлен в виде дерева попарных сложений.
В этой схеме попарно складываются наиболее близкие друг к другу числа, что влечет за собой минимизацию погрешности.
Интересный трюк использован в строках 12-18. Вычисление g теперь происходит в два действия КУС, и получаемый итоговый результат точнее, чем при одинарном умножении, эквивалентном с точки зрения школьных правил упрощения выражений. Этот приём я подсмотрел в открытой библиотеке Intel для сопроцессора 8087, опубликованной около 2000 года.
Прогоним тест и увидим, что погрешность теперь в пределах допустимой.
Другой хороший способ уменьшения погрешности – вычисление суммы с конца. Тогда и сумма, и очередное слагаемое растут одновременно, сводя к минимуму ошибку округления. Ещё немного усовершенствовав эту идею, придем к известной формуле Горнера.
Схема Горнера позволяет вычислить значение многочлена за минимальное количество арифметических операций, особенно если использовать КУС.
Получился очень компактный код, но мы утратили возможность внеочередного выполнения команд: каждое действие зависит от результата предыдущего.
Попробуем немного «испортить» код, вернув в него параллелизм. Сгруппируем слагаемые попарно, вынося за скобки общие множители.
Это похоже на схему Горнера, но не относительно x, а относительно x2.
Повторим группировку, теперь относительно x4 .
Получившаяся схема Эстрина требует ненамного больше действий, чем схема Горнера, при этом предоставляя возможность параллельного вычисления подвыражений.
Прогоним тест производительности с помощью Google Benchmark.
Дерево попарных сложений оказывается не только точнее, но и быстрее прямолинейного школьного подхода.
Схема Горнера, несмотря на наименьшее количество действий, оказалась самой медленной. Могу рекомендовать ее только для реализации в маломощных, например, в носимых устройствах, где количество операций, а значит, энергозатраты, в приоритете. Тем более, что из энергоэффективных устройств может быть вырезан блок внеочередного исполнения команд для экономии и места на кристалле, и энергии.
Ну а схема Эстрина, при своей неочевидности и непопулярности, рулит.
Ищите неочевидные подходы.