(5 оценок, среднее: 4,80 из 5)
Загрузка...
Учебная работа. Реферат: Варианты алгоритма возведения в степень повышение точности и ускорение
Максим М. Гумеров
Как, никто этого еще не выдумал?
Не берусь судить. возможно, задачка о том, как очень стремительно возвести действительное положительное число в произвольную действительную степень, решалась приблизительно настолько же нередко, как и вставала, — а вставала, полагаю, не раз. И все таки не так издавна я с страхом нашел, что RTL из состава Borland Delphi крайних версий (как Delphi 6, так и Delphi 7) подступает к решению не наиболее мастерски, чем прилежный пятиклассник, в первый раз столкнувшийся с таковой неувязкой.
Взглянем на начальный код функции Power из модуля Math, разлюбезно предоставленный Borland Software:
function Power(const Base, Exponent: Extended): Extended;
begin
if Exponent = 0.0 then
Result := 1.0 { n**0 = 1 }
else if (Base = 0.0) and (Exponent > 0.0) then
Result := 0.0 { 0**n = 0, n > 0 }
else if (Frac(Exponent) = 0.0) and (Абс(Exponent) <= MaxInt) then
Result := IntPower(Base, Integer(Trunc(Exponent)))
else
Result := Exp(Exponent * Ln(Base))
end;
Броско, что в благих целях оптимизации машина — комплекс технических средств, предназначенных для автоматической обработки информации в процессе решения вычислительных и информационных задач) (либо вычислительной системы) которое делает арифметические и логические операции данные программкой преобразования инфы управляет вычислительным действием и коор оставляют наедине с целой массой ветвлений, приводящих его, в конце концов, в общем случае к несчастному решению пятиклассника, а конкретно, к очевидной формуле
(1) x**y = exp(ln(x**y)) = exp(y*ln(x)).
тут x**y значит возведение x в степень y, a exp(x) = e**x.
Что отвратительного в таком подходе к решению? Во-1-х, в набор инструкций FPU не заходит ни операция вычисления exp(x), ни взятия натурального логарифма ln(x). Соответственно, итог рассчитывается в несколько шагов, в то время как можно пойти наиболее прямым методом, как будет показано ниже. Из-за этого падает скорость вычисления; не считая того, тут действует интуитивное правило, которое грубо можно сконструировать так: чем больше операций производится над числом с плавающей запятой в регистрах сопроцессора, тем больше будет и суммарная погрешность результата.
ПРИМЕЧАНИЕ
Позднейшая проверка показала, что как Visual C из Visual Studio .NET, так и C++ Builder 4.5 реализуют возведение в степень наиболее отменно. Применяемый в их вариант философски не различается от того решения, которое я желаю предложить.
Есть предложение
Давайте проведем инвентаризацию. Какие аннотации сопроцессора соединены с возведением в степень либо взятием логарифма? Приведу короткую выдержку из [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) значит вычисление квадратного корня.
Вот, в общем-то, и все. Но уже на 1-ый взор этого хватает, чтоб осознать, что задачка быть может решена наиболее прямо, чем дает 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) естественно выражаются таковым ассемблерным кодом:
;Во-1-х, вычисляем z=y*log2(x):
fld y ;Загружаем основание и показатель степени
fld x
fyl2x ;Стек FPU сейчас содержит: ST(0)=z
;Сейчас считаем 2**z:
fld st(0) ;Создаем еще одну копию z
frndint ;Округляем
fsubr st(0),st(1) ;ST(1)=z, ST(0)=z-trunc(z)
f2xm1 ;ST(1)=z, ST(0)=2**(z-trunc(z))-1
fld1
faddp ;ST(1)=z, ST(0)=2**(z-trunc(z))
fscale ;ST(1)=z, ST(0)=(2**trunc(z))*(2**(z-trunc(z)))=2**t
fxch st(1)
fstp st ;Итог остается на верхушке стека ST(0)
ПРЕДУПРЕЖДЕНИЕ
Перед выполнением этого фрагмента кода необходимо убедиться, что биты управления округлением в слове управления FPU установлены в режим округления к нулю. В Delphi это проще всего создать с помощью функции SetRoundMode (модуль Math):
SetRoundMode(rmTruncate);
ПРИМЕЧАНИЕ
Потому что на микропроцессорах Intel Pentium IV последовательное неоднократное переключение меж 2-мя (но не наиболее) состояниями слова управления FPU производится еще резвее, чем на ранешних моделях, можно рассчитывать, что даже в тех ситуациях, когда придется перемежать вызов этого фрагмента кода с действиями, требующими другого режима округления, при работе на современной технике это не приведет к чрезмерным доп временным затратам. Подробности см., к примеру, в [3].
Полный код работоспособной функции на Object Pascal смотрится так:
function _Power(const x,y:FLOATTYPE):FLOATTYPE; //x>0!
asm
fld y
fld x
fyl2x
fld st(0)
frndint
fsubr st(0),st(1)
f2xm1
fld1
faddp
fscale
fxch st(1)
fstp st
end;
СОВЕТ
Имеет смысл сделать перегруженные версии функции для разных типов аргументов FLOATTYPE, потому что на практике нередко основным недочетом интегрированной функции будет то, что она (как и все вызываемые ею функции) воспринимает в качестве аргументов действительные числа типа Extended, что приводит к очень значимым затратам на конвертирование форматов при загрузке в стек FPU.
Что мы достигнули?
Опыты проявили, что предложенный вариант функции возведения в степень увеличивает точность вычислений на один-два знака опосля запятой. Потому что создателю было несколько лень писать неспешный код для сверхточного возведения в степень с целью проверки точности предложенного метода, то количество же этих операций в предложенном методе и их общее число в реализации функций 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, — аппроксимация сплайнами. Говоря обычным языком (может быть, очень обычным – пусть меня извинят спецы), сплайн – это кривая, моделирующая форму, принимаемую упругим стержнем, деформированным методом закрепления в данных точках. Она проходит буквально через данные точки, подчиняясь при всем этом неким доп условиям, а именно, условию непрерывности 2-ой производной. Есть разные виды сплайнов. В данной для нас работе довольно удобно внедрение кубических сплайнов. Кубический сплайн на любом отрезке меж 2-мя поочередными (в порядке возрастания координаты 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 значений на отрицательной полуоси необходимы были лишь для того, чтоб отразить б (если для представления коэффициентов употребляется тип 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,[eax*4] //Другими словами eax:=i*4
add eax,eax // *2
add eax,1 // +1 = i*8+1 (дальше при доступе к массиву употребляется eax*4, другими словами i*32+4,
// т.к. любая строчка по 4*8=32 б и заполнитель сначала – 4 б.
// Если б не было заполнителя, последнюю аннотацию необходимо было бы убрать.
fadd st,st
fld1
fdivrp //=0.5
fmulp //st(1)=x, st(0)=0.5*trunc(2x)
fsubp //st(0)=dx
//Подсчет по схеме Горнера. Мне чудилось, что можно создать это резвее,
//пустив параллельно несколько цепочек вычислений, но пока это не удалось создать.
fld qword ptr coeffspow[4*eax]
fmul st,st(1)
fld qword ptr coeffspow[4*eax+8]
faddp
fmul st,st(1)
fld qword ptr coeffspow[4*eax+16]
faddp
fmul st,st(1)
fld qword ptr coeffspow[4*eax+24]
faddp
fxch st(1)
fstp st //Освобождаем ненадобный регистр
end;
ПРЕДУПРЕЖДЕНИЕ
Выполнение этого фрагмента изменяет содержимое регистра EAX.
Оценим погрешность приближения. Потому что итог, получаемый как _Power(2,x) (функция _Power приведена сначала статьи), заранее поточнее, чем Exp2(x), то в качестве оченки примем относительное отклонение значения крайней функции от значения первой: Epsilon=абс( Exp2(x) — _Power(2,x) ) / _Power(2,x). Очевидно, выражение имеет смысл, если _Power(2,x)<>0.
Если выстроить график относительной погрешности, становится видно, что в границах всякого из 1998 отрезков он имеет форму кривой с одним максимумом, сходящей к нулю на концах отрезка. При всем этом пределы колебаний величины погрешности остаются неизменными на всех отрезках, не считая нескольких крайних – на их погрешность растет. Если не принимать во внимание эти отрезки, и ограничить область допустимых значений аргумента числом 990 (т.е. x<990), то для описания поведения относительной погрешности зависимо от x довольно показать ее график на 2-ух крайних допустимых для значений 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; //абс(y*log2(x))<990
asm
fld y
fld x
fldz //Сравним основание степени
fcomip st,st(1) // с 0 и соответственно установим флаги микропроцессора
je @Zero
FYL2X //Стек: ST(0)=t=y*log2(x)
fldz
fcomip st,st(1) //Флаги выставляем соответственно числу 0-y*log2(x)
ja @Reverse //Если 0>y*log2(x), то сосчитаем 2**|y*log2(x)|, а опосля инвертируем
call Core_Exp2
jmp @Exit
@Zero:
fxch st(1)
fstp st //Освобождаем ненадобный регистр
jmp @Exit
@Reverse:
fabs //Берем абсолютную величин
call Core_Exp2
fld1 //Считаем оборотное
fdivrp //1/(2**|y*log2(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; // абс(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%) даже при весьма огромных показателях:
Набросок 5.
Заключение
Таковым образом, нам удалось получить функции, превосходящие встроенную по скорости от 2-ух до 4 раз при погрешности порядка 0.004% — 0.01%. Не исключено, что существует возможность провести наиболее доброкачественную и наиболее прибыльную в смысле временных издержек аппроксимацию функций; может быть, даже по другому принципу, а не так, как это было изготовлено тут (т.е. исходя из соотношения x**y=2**(y*log2(x)) ).
Для тех же случаев, когда нужна высочайшая точность вычислений, в качестве первого камня фундамента была рассмотрена функция, исправляющая недочет Delphi RTL. Непременно, это направление также достойно последующих исследовательских работ с целью убыстрить заранее неспешные вычисления с завышенной точностью.
Перечень
литературы
Intel® Architecture Software Developer’s Manual: том 2, Instruction set reference. Можно отыскать на веб-сайте Intel www.intel.com.
Intel® VTune™ Performance Analyzer, гипертекстовая справка. И совершенно, VTune – превосходный инструмент для поиска шероховатостей в программке.
Intel® Pentium® 4 and Intel® Xeon™ Processor Optimization Reference Manual. Все там же, на веб-сайте Intel.
]]>