Варианты алгоритма возведения в степень повышение точности и ускорение
Варианты алгоритма возведения в степень: повышение точности и ускорение.
Максим М. Гумеров
Как, никто этого еще не придумал?
Не берусь судить. Вероятно, задача о том, как максимально быстро возвести действительное положительное число в произвольную действительную степень, решалась примерно столь же часто, как и вставала, - а вставала, полагаю, не раз. И все же не так давно я с ужасом обнаружил, что 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)
Наверняка у вас есть товары или услуги, продажа которых приносит вам максимальную прибыль. Для быстрого старта в сети вам необходимо создание посадочной страницы (одностраничного сайта), на которой будет размещена информация о маржинальных товарах/услугах интернет магазина. За 8 лет опыта разработки конверсионных страниц мы выработали оптимальную структуру, которая позволит привлекать через landing page больше продаж. На такую структуру «одевается» ваш контент — фирменный стиль, тексты, фотографии, уникальные торговые предложения, после чего страница выходит в свет. Разработка лендинга и запуск в сети — до 7 рабочих дней. Стоит отметить, что в разработку самой посадочной страницы входит и написание копирайтером продающих текстов для вашего бизнеса, чтобы каждый посетитель страницы захотел совершить покупку именно у вас. Результат: качественно разработаная продающая посадочная страница, которая готова приносить вам новых клиентов. © 2016 - 2022 BigEdu.ru
|