Варианты алгоритма возведения в степень повышение точности и ускорение
Варианты алгоритма возведения в степень: повышение точности и ускорение.
Максим М. Гумеров
Как, никто этого еще не придумал?
Не берусь судить. Вероятно, задача о том, как максимально быстро возвести действительное положительное число в произвольную действительную степень, решалась примерно столь же часто, как и вставала, - а вставала, полагаю, не раз. И все же не так давно я с ужасом обнаружил, что 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 (Abs(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.
Что плохого в таком подходе к решению? Во-первых, в набор инструкций 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) означает вычисление квадратного корня.
Вот, в общем-то, и все. Но уже на первый взгляд этого хватает, чтобы понять, что задача может быть решена более прямо, чем предлагает 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) естественно выражаются таким ассемблерным кодом:
| ;Во-первых, вычисляем 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 последовательное многократное переключение между двумя (но не более) состояниями слова управления FPU выполняется гораздо быстрее, чем на ранних моделях, можно рассчитывать, что даже в тех ситуациях, когда придется перемежать вызов этого фрагмента кода с действиями,
Оценок: 1005 (Средняя 5 из 5)
Одними из наиболее популярных услуг на рынке IT-технологий являются создание и продвижение лендингов. Они способны положительно влиять на деятельность любого бизнес-проекта в интернете. Судя по многочисленным отзывам, заказавшие создание лендингов люди ни разу не пожалели о потраченных деньгах. Они вложили в будущее, которое неразрывно связано с интернетом. Всё больше и больше предпринимателей обращаются к услугам разных агентств, веб-студий, чтобы заказать создание лендинга у профессионалов. © 2017 - 2022 ReferatWorld.ru
|