Варианты алгоритма возведения в степень повышение точности и ускорение (стр. 2 из 3)

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

array [0..1998] of packed record c3,c2,c1,c0:double end;

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

//Предполагаю, что выставлена опция компилятора {$Align 8} Type TArr=packed record Padding:integer; //Фиктивный 4-байтовый заполнитель - чтобы массив выравнялся по 8 байтам C:array [0..1998] of packed record c3,c2,c1,c0:double end; //Собственно коэффициенты end;

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

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

ПРЕДУПРЕЖДЕНИЕ Как и для предыдущей функции, нужно обеспечить установку бит управления округлением в режим округления к нулю.

function Exp2(x:FLOATTYPE):FLOATTYPE; //0<=x<999 asm fld x call Core_Exp2 //Оформим основную часть в виде процедуры, т.к. она будет использоваться не только здесь - // - да и перегрузку функций для другого типа аргумента так делать удобнее. end; procedure Core_Exp2; //На вершине стека FPU находится аргумент var i:integer; //Сюда получим индекс в массиве asm fld st //Копируем аргумент fadd st,st //st(1)=x, st(0)=2x fistp i //Достаем i (индекс равен trunc(2x)); st(0)=x fild i //Полагаемся на т.н. Store-Forwarding: округленное значение передается сразу инструкции // fild, не ожидая, пока данные будут записаны в память; st(1)=x, st(0)=trunc(2x) mov eax,i fld1 //st(2)=x, st(1)=trunc(2x), st(0)=1 lea eax,[eax4] //То есть eax:=i4 add eax,eax // 2 add eax,1 // +1 = i8+1 (далее при доступе к массиву используется eax4, то есть i32+4, // т.к. каждая строка по 48=32 байта и заполнитель в начале – 4 байта. // Если бы не было заполнителя, последнюю инструкцию нужно было бы убрать. fadd st,st fld1 fdivrp //=0.5 fmulp //st(1)=x, st(0)=0.5trunc(2x) fsubp //st(0)=dx //Подсчет по схеме Горнера. Мне казалось, что можно сделать это быстрее, //пустив параллельно несколько цепочек вычислений, но пока это не удалось сделать. fld qword ptr coeffspow[4eax] fmul st,st(1) fld qword ptr coeffspow[4eax+8] faddp fmul st,st(1) fld qword ptr coeffspow[4eax+16] faddp fmul st,st(1) fld qword ptr coeffspow[4eax+24] faddp fxch st(1) fstp st //Освобождаем ненужный регистр end;
ПРЕДУПРЕЖДЕНИЕ Выполнение этого фрагмента изменяет содержимое регистра EAX.

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

Если построить график относительной погрешности, становится видно, что в пределах каждого из 1998 отрезков он имеет форму кривой с одним максимумом, сходящей к нулю на концах отрезка. При этом пределы колебаний величины погрешности остаются постоянными на всех отрезках, кроме нескольких последних – на них погрешность возрастает. Если не принимать во внимание эти отрезки, и ограничить область допустимых значений аргумента числом 990 (т.е. x<990), то для описания поведения относительной погрешности в зависимости от x достаточно показать ее график на двух последних допустимых для значений x отрезках:

Рисунок 1. Максимальная погрешность приближения функции Exp2=2**x (при x менее 990) не превышает 0,004%.

СОВЕТ Мы отсекли отрезки, лежащие правее точки x=990. Следовательно, размер таблицы коэффициентов можно несколько сократить: индекс последнего элемента должен быть 990*2=1980, а не 1998. “Лишние” 19 последних строк таблицы можно просто удалить. Логично также изменить текст комментария в начале функции Exp2.

Новый вариант функции возведения в степень

Изменим реализацию возведения в степень в соответствии с предложенной аппроксимацией для 2**x:

function New_Power(x,y:FLOATTYPE):FLOATTYPE; //abs(ylog2(x))<990 asm fld y fld x fldz //Сравним основание степени fcomip st,st(1) // с 0 и соответственно установим флаги процессора je @Zero FYL2X //Стек: ST(0)=t=ylog2(x) fldz fcomip st,st(1) //Флаги выставляем соответственно числу 0-ylog2(x) ja @Reverse //Если 0>ylog2(x), то сосчитаем 2*\|ylog2(x)\|, а после инвертируем call Core_Exp2 jmp @Exit @Zero: fxch st(1) fstp st //Освобождаем ненужный регистр jmp @Exit @Reverse: fabs //Берем абсолютную величин call Core_Exp2 fld1 //Считаем обратное значение: fdivrp //1/(2*\|ylog2(x)\|) @Exit: end;
ПРЕДУПРЕЖДЕНИЕ В этом фрагменте используется инструкция FCOMIP, впервые появившаяся на процессорах Pentium Pro. Любителям антиквариата придется использовать вместо пары команд FCOMIP / JE блок

FCOMP FSTSW TEST AX, 16384 JNZ @Zero //Вместо je @Zero
ПРЕДУПРЕЖДЕНИЕ А вместо FCOMIP / JA - блок

FCOMP FSTSW TEST AX, 256 or 16384 //0<= y*log2(x) ? JZ @Reverse //Нет, случай со взятием обратного значения
ПРЕДУПРЕЖДЕНИЕ Вдобавок в этом случае изменяется регистр EAX.

Результаты тестирования отражены на графиках:

Рисунок 2. Временные затраты: New_Power – новая функция, Power – из состава RTL Borland Delphi.

Подпись X-0.511 на оси абсцисс отражает тот факт, что при проведении испытаний брались значения целые значения X, к которым затем прибавлялось число 0.511, чтобы гарантировать, что основание степени – число нецелое (т.е. чтобы рассматривать по возможности общий случай).

Черная линия поверх красного набора – сглаженные временные затраты функции Power, фиолетовая поверх синего – функции New_Power.

Замеры временных затрат производились с помощью инструкции RDTSC (процессоры начиная с Pentium):

function time:int64; //Вспомогательная функция для подсчета времени работы asm rdtsc end;

и далее в коде

t:=time(); ... writeln(time()-t);

RDTSC возвращает в регистровой паре EDX:EAX число тактов процессора, прошедших с момента последнего сброса (reset). Машинный код инструкции – 0Fh, 31h.

Относительная погрешность ведет себя на удивление стабильно, изменяясь в пределах от 0 до 0,0040%. Поэтому достаточно представительным множеством значений аргумента является, к примеру, промежуток (0, 1000).

Рисунок 3.

Видно, что оцененная относительная погрешность (фактически - отклонение от значения, возвращаемого встроенной функцией) на самом деле не превосходит 0.004% !

В случае показателя степени 17 колебания становятся намного чаще, однако общая картина та же.

Аппроксимация функции log2x и “специализация” возведения в степень

Логарифмирование плохо поддается аппроксимации с помощью кубических сплайнов – точнее, мне удалось это сделать, причем с весьма высокой точностью, но лишь ценой проигрыша по времени в сравнении с использованием FYL2X. Однако здесь есть что предпринять и не прибегая к сплайнам.

Как известно, функция ln(1+x) при |x|<1 разлагается в ряд Тейлора следующим образом:

ln(1+x)=x-x2/(1*2)+x3/(1*2*3)+…+ xi/i!+…

Если абсолютная величина x достаточно мала, члены ряда, уже начиная с третьего, достаточно слабо сказываются на результате. Поэтому для значений x, достаточно близких к 1 (чтобы остаться в оговоренных выше рамках приемлемых погрешностей, x должен отстоять от 1 не больше чем на 0.01), вычисление log2(x)=ln(x)/ln(2)=ln(x)*log2(e)=ln(1+(x-1))*log2(e) можно заменить вычислением (t-t*t/2)*log2(e), где t=x-1.

Это позволяет построить еще один вариант функции возведения в степень для значений основания, близких к 1. В нем нет инструкции FYL2X, а вместо нее присутствует блок инструкций, помеченных символом “ * ” (знак “~” означает приближенное равенство):

function New_Power_XNear1(x,y:FLOATTYPE):FLOATTYPE; // abs(y*log2(x))<990 asm fld y fld x fldz fcomip st,st(1) je @Zero fld1 (*) fsub st(1),st (*) fld st(1) (*) //st(0)=1; st(1)=st(3)=t=x-1, st(2)=1, st(4)=y fld1 (*) fadd st,st (*) fdivp st(2),st (*) //st(0)=st(2)=t, st(1)=1/2, st(3)=y fmul st,st (*) fmulp st(1),st (*) //st(0)=1/2*t*t, st(1)=t, st(2)=y fsubp st(1),st (*) //st(0)=t-t*t/2 ~ ln(x), st(1)=y fldl2e (*) //Загружаем константу log2(e) fmulp (*) //st(0)~log2(x), st(1)=y fmulp (*) //st(0)~y*log2(x) fldz fcomip st,st(1) ja @Reverse call Core_Exp2 jmp @Exit @Zero: fxch st(1) fstp st //Освобождаем ненужный регистр jmp @Exit @Reverse: fabs call Core_Exp2 fld1 fdivrp @Exit: end;

Таким образом, нам в этом случае (при x, близких к 1) удается избавиться от всех инструкций FPU, принадлежащих к группе трансцендентных, что приводит к впечатляющему росту производительности:

Рисунок 4. Временные затраты: New_Power_XNear1 – специализированный вариант New_Power.

К сожалению, с ростом показателя степени максимальная погрешность растет, оставаясь, впрочем, в оговоренных пределах (т.е. меньше 0,1%; более того – меньше 0,01%) даже при очень больших показателях: