Варианты алгоритма возведения в степень: повышение точности и ускорение

Максим М. Гумеров 

Как, никто этого еще не придумал?

Не берусь судить. Вероятно, задача о том, как максимально быстро возвести действительное положительное число в произвольную действительную степень, решалась примерно столь же часто, как и вставала, - а вставала, полагаю, не раз. И все же не так давно я с ужасом обнаружил, что RTL из состава Borland Delphi последних версий (как Delphi 6, так и Delphi 7) подходит к решению не более профессионально, чем прилежный пятиклассник, впервые столкнувшийся с такой проблемой.

Взглянем на исходный код функции Power из модуля Math, любезно предоставленный Borland Software:

Примечательно, что в благих целях оптимизации процессор оставляют наедине с целой толпой ветвлений, приводящих его, в конце концов, в общем случае к пресловутому решению пятиклассника, а именно, к тривиальной формуле

(1) x**y = exp(ln(x**y)) = exp(y*ln(x)).

Здесь x**y означает возведение x в степень y, a exp(x) = e**x.

Что плохого в таком подходе к решению? Во-первых, в набор инструкций FPU не входит ни операция вычисления exp(x), ни взятия натурального логарифма ln(x). Соответственно, результат вычисляется в несколько этапов, в то время как можно пойти более прямым путем, как будет показано ниже. За счет этого падает скорость вычисления; кроме того, здесь действует интуитивное правило, которое грубо можно сформулировать так: чем больше операций выполняется над числом с плавающей запятой в регистрах сопроцессора, тем больше будет и суммарная погрешность результата.

Есть предложение

Давайте проведем инвентаризацию. Какие инструкции сопроцессора связаны с возведением в степень или взятием логарифма? Приведу краткую выдержку из [1] и [2]:

F2XM1 – вычисляет 2**x-1, где -1<x<1.

FSCALE (масштабирование степенью двойки) - фактически считает 2**trunc(x), где trunc(x) означает округление к 0, т.е. положительные числа округляются в меньшую сторону, отрицательные – в большую.

FXTRACT – извлекает мантиссу и экспоненту действительного числа.

FYL2X – вычисляет y*log2(x), где log2(x) – логарифм x по основанию 2.

FYL2XP1 – вычисляет y*log2(x+1) для -(1-1/sqrt(2))<x<1-1/sqrt(2) c большей точностью, нежели FYL2X. Здесь sqrt(x) означает вычисление квадратного корня.

Вот, в общем-то, и все. Но уже на первый взгляд этого хватает, чтобы понять, что задача может быть решена более прямо, чем предлагает RTL Borland Delphi.

Действительно, почему не заменить показатель степени в (1) на 2? Ведь неперово число отнюдь не является родным для двоичной арифметики! Тогда получится

(2) x**y = 2**log2(x**y) = 2**(y*log2(x)).

Это выражение для x**y в соответствии с вышеозначенными пятью инструкциями можно представить как композицию функций в таком виде:

(3) f(z)=2**z,

(4) g(x,y)=y*log2(x),

(5) xy =f(g(x,y)).

Так как вычислить f(z) в одно действие невозможно, приходится считать так:

(6) f(z)=2**z=2**(trunc(z)+(z-trunc(z)))=(2**trunc(z)) * (2**(z-trunc(z))).

Формулы (4)-(6) естественно выражаются таким ассемблерным кодом:

Полный код работоспособной функции на Object Pascal выглядит так:

Чего мы достигли?

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

К сожалению, выигрыш в скорости абсолютно не ощущается. Это вполне объяснимо: согласно приложению C (“IA-32 Instruction Latency and Throughput”) документа [3], из всего этого фрагмента основная вычислительная нагрузка ложится на трансцендентные (ответственность за не вполне корректное применение термина ложится не на меня, а на господ из Intel) операции, а именно – FYL2X, FRNDINT, F2XM1 и FSCALE. Количество же этих операций в предложенном алгоритме и их общее число в реализации функций ln(x) и exp(x) в RTL Delphi одинаково.

Конечно, хотелось бы увеличить и скорость вычислений. Но мир не идеален, и за это придется расплачиваться все той же точностью. Как правило, в каждой ситуации существует предел допустимых погрешностей при расчетах. В иллюстративных целях я задался максимальной допустимой относительной погрешностью 0,0001=0,1%. В действительности, как будет видно из графиков относительной погрешности, удалось достичь еще большей точности.

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

Аппроксимация функции 2x

Эта мера позволит нам избавиться сразу и от длительной F2XM1, и от выполняющейся ненамного быстрее FSCALE.

Существует бесконечное множество способов приблизить функцию f(x). Один из наиболее простых в вычислительном плане – подбор подходящего по точности многочлена g(x)=anxn+an-1xn-1+...+a1x+a0. Его коэффициенты могут быть постоянны, а могут некоторым образом зависеть от x. В первом случае коэффициенты легко найти методом наименьших квадратов, взяв значения исходной функции в нескольких точках и подобрав коэффициенты так, чтобы в этих точках многочлен принимал значения, как можно более близкие к значениям функции (подробное описание полиномиальной аппроксимации функций и метода наименьших квадратов можно найти в книгах, посвященных курсам вычислительной математики или обработке экспериментальных данных). Простота метода оборачивается существенным недостатком: он подчас неплох для выявления качественных тенденций, но плохо отражает исходную функцию количественно, причем, как правило, погрешность растет с уменьшением степени многочлена n, а скорость вычисления g(x) с ростом n падает.

Достойная альтернатива, позволяющая достаточно точно приближать гладкие кривые, такие, как y=2**x, - аппроксимация сплайнами. Говоря простым языком (возможно, чересчур простым – пусть меня извинят специалисты), сплайн – это кривая, моделирующая форму, принимаемую упругим стержнем, деформированным путем закрепления в заданных точках. Она проходит точно через заданные точки, подчиняясь при этом некоторым дополнительным условиям, в частности, условию непрерывности второй производной. Существуют различные виды сплайнов. В этой работе достаточно практично использование кубических сплайнов. Кубический сплайн на каждом отрезке между двумя последовательными (в порядке возрастания координаты x) эталонными точками (x,f(x)) описывается полиномом третьей степени g(x)=a3x3+a2x2+a1x+a0, где набор коэффициентов (a0,a1,a2,a3) свой для каждого отрезка. Поиск этих коэффициентов – не слишком сложная задача, но описание метода ее решения выходит за рамки этой статьи. Таблица коэффициентов, получающаяся после учета всех замечаний этого раздела, прилагается к статье.

Итак, я ограничусь лишь использованием полученных мною значений коэффициентов. Чтобы обеспечить необходимую точность на промежутке 0<=x<999, мне понадобились в качестве эталонных 2039 точек, которым соответствовали значения x=(i-40)/2, i=0..2038. Сорок значений на отрицательной полуоси нужны были только для того, чтобы отразить поведение сплайна в этой части плоскости, слегка скорректировав таким образом его вид на остальных отрезках; в вычислениях эти 40 отрезков не участвуют, т.к. для значений x<0 можно воспользоваться (без ощутимого проигрыша в скорости или точности) соотношением 2**(-|x|)=1/(2**|x|).

Итак, у нас есть таблица коэффициентов, представленная в виде массива из 1999 блоков по 8*4 байт (если для представления коэффициентов используется тип double). На Object Pascal такой массив описывается типом

На практике возникает тонкий момент. Дело в том, что Delphi почему-то отказывается выравнивать массивы Double’ов по границе 8 байт. Лично у меня получается так, что адрес первого элемента всегда бывает кратен 4, но никогда – 8. Поэтому перед началом массива я вставляю заполнитель, чтобы избежать медленного чтения некоторых double’ов, которые частично лежат в одной 64- или 32-байтной линейке кэша, а частично – в следующей:

На вход функции Exp2 поступает единственный аргумент x - возводимое в степень число. Как можно реализовать функцию?

Вот как это сделал я.

Оценим погрешность приближения. Так как результат, получаемый как _Power(2,x) (функция _Power приведена в начале статьи), заведомо точнее, чем Exp2(x), то в качестве оченки примем относительное отклонение значения последней функции от значения первой: Epsilon=abs( Exp2(x) - _Power(2,x) ) / _Power(2,x). Разумеется, выражение имеет смысл, если _Power(2,x)<>0.

Если построить график относительной погрешности, становится видно, что