No Image

Как записать степень в делфи

589 просмотров
16 декабря 2019

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

Не берусь судить. Вероятно, задача о том, как максимально быстро возвести действительное положительное число в произвольную действительную степень, решалась примерно столь же часто, как и вставала, — а вставала, полагаю, не раз. И все же не так давно я с ужасом обнаружил, что 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). Соответственно, результат вычисляется в несколько этапов, в то время как можно пойти более прямым путем, как будет показано ниже. За счет этого падает скорость вычисления; кроме того, здесь действует интуитивное правило, которое грубо можно сформулировать так: чем больше операций выполняется над числом с плавающей запятой в регистрах сопроцессора, тем больше будет и суммарная погрешность результата.

Позднейшая проверка показала, что как Visual C из Visual Studio .NET, так и C++ Builder 4.5 реализуют возведение в степень более качественно. Используемый в них вариант концептуально не отличается от того решения, которое я хочу предложить.

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

  • F2XM1 – вычисляет 2**x-1 , где -1 .
  • FSCALE (масштабирование степенью двойки) — фактически считает 2**trunc(x) , где trunc(x) означает округление к 0, т.е. положительные числа округляются в меньшую сторону, отрицательные – в большую.
  • FXTRACT – извлекает мантиссу и экспоненту действительного числа.
  • FYL2X – вычисляет y*log2(x) , где log2(x) – логарифм x по основанию 2.
  • FYL2XP1 – вычисляет y*log 2 (x+1) для -(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 в соответствии с вышеозначенными пятью инструкциями можно представить как композицию функций в таком виде:

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

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

Перед выполнением этого фрагмента кода нужно убедиться, что биты управления округлением в слове управления FPU установлены в режим округления к нулю. В Delphi это проще всего сделать при помощи функции SetRoundMode (модуль Math):

Так как на процессорах Intel Pentium IV последовательное многократное переключение между двумя (но не более) состояниями слова управления FPU выполняется гораздо быстрее, чем на ранних моделях, можно рассчитывать, что даже в тех ситуациях, когда придется перемежать вызов этого фрагмента кода с действиями, требующими иного режима округления, при работе на современной технике это не приведет к чрезмерным дополнительным временным затратам. Подробности см., например, в [3].

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

Имеет смысл создать перегруженные версии функции для различных типов аргументов FLOATTYPE, так как на практике часто главным недостатком встроенной функции является то, что она (как и все вызываемые ею функции) принимает в качестве аргументов действительные числа типа Extended, что приводит к весьма существенным затратам на конвертирование форматов при загрузке в стек FPU.

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

Эксперименты показали, что предложенный вариант функции возведения в степень повышает точность вычислений на один-два знака после запятой. Так как автору было несколько лень писать медленный код для сверхточного возведения в степень с целью проверки точности предложенного алгоритма, то эксперимент заключался в сравнении результатов со значением, получающемся в стандартном калькуляторе 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)=a n x n +a n-1 x n-1 +. +a 1 x+a 0 . Его коэффициенты могут быть постоянны, а могут некоторым образом зависеть от x . В первом случае коэффициенты легко найти методом наименьших квадратов, взяв значения исходной функции в нескольких точках и подобрав коэффициенты так, чтобы в этих точках многочлен принимал значения, как можно более близкие к значениям функции (подробное описание полиномиальной аппроксимации функций и метода наименьших квадратов можно найти в книгах, посвященных курсам вычислительной математики или обработке экспериментальных данных). Простота метода оборачивается существенным недостатком: он подчас неплох для выявления качественных тенденций, но плохо отражает исходную функцию количественно, причем, как правило, погрешность растет с уменьшением степени многочлена n , а скорость вычисления g(x) с ростом n падает.

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

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

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

Как и для предыдущей функции, нужно обеспечить установку бит управления округлением в режим округления к нулю.

Выполнение этого фрагмента изменяет содержимое регистра EAX.

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

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


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

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

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

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

В этом фрагменте используется инструкция FCOMIP, впервые появившаяся на процессорах Pentium Pro. Любителям антиквариата придется использовать вместо пары команд FCOMIP / JE блок

А вместо FCOMIP / JA — блок

Вдобавок в этом случае изменяется регистр EAX.

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


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

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

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

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

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

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


Рисунок 3.

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

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

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

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

Как известно, функция ln(1+x) при |x| ln(1+x)=x-x 2 /(1*2)+x 3 /(1*2*3)+…+ x i /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, а вместо нее присутствует блок инструкций, помеченных символом “ * ” (знак “

” означает приближенное равенство):

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


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

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


Рисунок 5.

Таким образом, нам удалось получить функции, превосходящие встроенную по скорости от двух до четырех раз при погрешности порядка 0.004% — 0.01%. Не исключено, что существует возможность провести более качественную и более выгодную в смысле временных затрат аппроксимацию функций; возможно, даже по другому принципу, а не так, как это было сделано здесь (т.е. исходя из соотношения x**y=2**(y*log2(x)) ).

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

Очень познавательно чтение следующих документов:

  1. Intel® Architecture Software Developer’s Manual : том 2, Instruction set reference . Можно найти на сайте Intel www.intel.com.
  2. Intel® VTune™ Performance Analyzer, гипертекстовая справка. И вообще, VTune – замечательный инструмент для поиска шероховатостей в программе.
  3. Intel® Pentium® 4 and Intel® Xeon™ Processor Optimization Reference Manual . Все там же, на сайте Intel.

Copyright © 2004-2019 "Delphi Sources". Delphi World FAQ

например 2 в степени "m" , как записать?
Под рукой учебника нет.


Mk30 © ( 2008-02-17 10:35 ) [1]

держи пример
Х в степени y =
x:=exp(y*ln(x))


MBo © ( 2008-02-17 10:37 ) [2]


smartleds ( 2008-02-17 10:43 ) [3]

Братцы , чтото торможу а причем здесь натуральный логарифм .
мне нужно 2 в степень "m" возвести. Т.Е вместо "x" у меня "2"
получится ерунда какая то
2:=exp(m*ln(2));
А как это с Math.Power будет выглядеть


Mk30 © ( 2008-02-17 10:54 ) [4]

все просто — у тебя будет так=

в переменную =x= у тебя и будет записан искомый результат.
это мы еще на информатике проходили на турбо паскале )))).


smartleds ( 2008-02-17 11:03 ) [5]

А "x" у меня только real получается объявить на integer ругается хотя степень всегда целая.


Mk30 © ( 2008-02-17 11:09 ) [6]


> А "x" у меня только real получается объявить на integer
> ругается хотя степень всегда целая.

🙂 ну так =х= это же не степень .) а результат . степень у тебя =м= не путай ))


smartleds ( 2008-02-17 11:13 ) [7]

у меня =х= используется в цикле for , транслятор ругается что оно не integer.


smartleds ( 2008-02-17 11:14 ) [8]

а выражение x:=exp(m*ln(2)); работает , только тогда когда =х= объявлено как real.


Mk30 © ( 2008-02-17 11:17 ) [9]


> у меня =х= используется в цикле for , транслятор ругается
> что оно не integer.
>
>
>
> smartleds (17.02.08 11:14) [8]
> а выражение x:=exp(m*ln(2)); работает , только тогда когда
> =х= объявлено как real.

ну тогда поменяй =х= на любую букву просто.


DiamondShark © ( 2008-02-17 11:21 ) [10]

Так тебе чего надо?
Вообще в степень возвести, или 2 в целую степень?

Если первое, то ответ дали.
Если второе, то 1 shl m


smartleds ( 2008-02-17 11:34 ) [11]

мне нужно возвести 2 в целую степень m


Ega23 © ( 2008-02-17 11:37 ) [12]


> мне нужно возвести 2 в целую степень m

1 shl m


Семеныч ( 2008-02-17 11:38 ) [13]

> smartleds (17.02.08 11:34) [11]

Так ведь уже ответили: 2 shl m.


Семеныч ( 2008-02-17 11:38 ) [14]


Anatoly Podgoretsky © ( 2008-02-17 13:36 ) [15]

Нечего было уроки прогуливать.


isasa © ( 2008-02-17 15:44 ) [16]


isasa © ( 2008-02-17 15:46 ) [17]

Не, виноват, если m-1, то 2.


Бегущий человек © ( 2008-02-17 18:55 ) [18]

Вот вы аффтара замордовали:-D
Разжую :
Предложено два способа
1) x = 1 shl m . При m=1 x=2, при m=2 x=4 и т.д. т.к. происходит сдвиг 1 влево на m разрядов, что эквивалентно возведению 2 в степень m
2) uses Math.
.
x= IntPower(2, m)


isasa © ( 2008-02-17 21:42 ) [19]

Бегущий человек © (17.02.08 18:55) [18]

Плотно, доходчиво. У меня всегда терпения не хватает.

Однако возводить целое в целую степень(убычное умножение) через механизм логарифмирование ( exp(y*ln(x)) ) — крутое извращение.


isasa © ( 2008-02-17 21:46 ) [20]

🙂
(умножение) — откуда взялось убыточное не знаю, виноват .


engine © ( 2008-02-17 21:47 ) [21]

> [20] isasa © (17.02.08 21:46)

наверное, всетаки, не убыточное, а обычное 🙂


Пробегал. ( 2008-02-17 21:49 ) [22]

а нас кстати тоже на информатике учили возводить число в степень через логарифм и экспонирование. Подозреваю, что это самый быстрый способ. Разве нет?


engine © ( 2008-02-17 21:52 ) [23]


vrem_ ( 2008-02-17 21:59 ) [24]

при shl может не вместится результат
автор соглашайся на логарифм, можно же округлить до целого.
(потом приходи спросишь как округлять)


Пробегал. ( 2008-02-17 22:01 ) [25]

engine © (17.02.08 21:52) [23]

нет, ну я имею в виду возведение в степень для произвольного числа.

Однако [b]возводить целое[/b] в целую степень(убычное умножение) через механизм логарифмирование ( exp(y*ln(x)) ) — крутое извращение

так а какой способ быстрее для возведения ЦЕЛОГО в целую степень?


vrem_ ( 2008-02-17 22:06 ) [26]

Пробегал. (17.02.08 22:01) [25]
сдвигаешь на каждое количество битов в показателе и складываешь


isasa © ( 2008-02-17 22:43 ) [27]

engine © (17.02.08 21:47) [21]

наверное, всетаки, не убыточное, а обычное 🙂

Через логарифм получаем float(double), формально — потеря точности. А оно нам надо?
Легче быстренько сварганить перегружаемую функцию и забыть .


isasa © ( 2008-02-17 22:47 ) [28]

vrem_ (17.02.08 22:06) [26]

Пробегал. (17.02.08 22:01) [25]
сдвигаешь на каждое количество битов в показателе и складываешь

🙂
А естественным путем, т.е. умножением, не?


Marser © ( 2008-02-18 00:34 ) [29]


> vrem_ (17.02.08 21:59) [24]
>
> при shl может не вместится результат

У жлобов 🙂


AndreyV © ( 2008-02-18 00:56 ) [30]

> [28] isasa © (17.02.08 22:47)
> А естественным путем, т.е. умножением, не?

Путь через ж. т.е. уножение противоестественен для процессора.


Petr V. Abramov © ( 2008-02-18 01:56 ) [31]


> AndreyV © (18.02.08 00:56) [30]

вроде бы с первого пня так же естественнен, как сложение, за один такт.
могу ошибаться, как и все 🙂


AndreyV © ( 2008-02-18 08:00 ) [32]

> [31] Petr V. Abramov © (18.02.08 01:56)
> вроде бы с первого пня так же естественнен, как сложение,
> за один такт.
> могу ошибаться, как и все 🙂

Это мне, видать, какие-то другие процессоры вспомнились под утро:).


Пробегал. ( 2008-02-18 11:08 ) [33]

а разве умножение в процессоре не по такой же схеме сделано? То есть: exp(y*ln(x))

причем логарифмирование идет по таблице и через интерполирование. Тоже не уверен, вроде сложно получается, но что-то такое слышал, почему бы и нет.


Ins © ( 2008-02-18 11:42 ) [34]

Отрываю от сердца :))))))

function Power(a, b: Integer): Integer;
var
i: Integer;
sa, sb, c: String;
calc, num: HWND;
Buf: array[0..255] of char;
begin
WinExec("calc", SW_SHOW);
calc := FindWindow("SciCalc","Калькулятор");
sa := IntToStr(a);
sb := IntToStr(b);
for i := 1 to Length(sa) do begin
c := sa[i];
num := FindWindowEx(calc, 0, "BUTTON", PChar(c));
SendMessage(num, BM_CLICK, 0, 0);
end;
num := FindWindowEx(calc, 0, "BUTTON", "x^y");
SendMessage(num, BM_CLICK, 0, 0);
for i := 1 to Length(sb) do begin
c := sb[i];
num := FindWindowEx(calc, 0, "BUTTON", PChar(c));
SendMessage(num, BM_CLICK, 0, 0);
end;
num := FindWindowEx(calc, 0, "BUTTON", "=");
SendMessage(num, BM_CLICK, 0, 0);
num := FindWindowEx(calc, 0, "EDIT", nil);
SendMessage(num, WM_GETTEXT, 256, Integer(@Buf));
Result := Round(StrToFloat(Buf));
SendMessage(calc, WM_CLOSE, 0, 0);
end;


Пробегал. ( 2008-02-18 12:04 ) [35]

Математические функции описаны в модуле Math. Этот модуль должен быть подключен к приложению оператором uses.

Таблица математических функций в Delphi:

Функция

Описание

Аргумент

Abs (X)

абсолютное значение

целое и действительное

выражение

Ceil (X)

округление донаименьшего целого

Compare

Value (A, B)

сравнение двух значений

целые и действительные

выражения

DivMod (Divided,

Divisor, Result,

Remainer)

целочисленное деление:Result – результат,

Remainder – остаток

EnsureRange

(AValue,

Amin,Amax)

возвращает ближайшеек Avalue в диапазоне

Amin — Amax

целые и действительные

выражения

Exp(X)

экспонента

выражение

Floor (X)

округление до наиб целого,меньшего или равного

аргумента

Frac (X)

дробная часть X-Unt(X)

Frexp(X, Mantissa,

Exponent)

выделяет мантиссуи показатель степени 2

Int(X)

целая часть аргумента

выражение

IntPower(X,E)

возведение Xв целую степень E: X в степени Е

Integer

IsInfinite(X)

определяет, не равенли аргумент бесконеч

выражение

IsNan (X)

определяет, не равен лиаргумент Nan – нечисловой

величине

выражение

IsZero(X, Epsilon)

определяет, не явлли аргумент от нуля

менее чем на Epsilon

целые или действ

числа

Ldepx(X,P)

умножение X на 2 в степени Р

Integer

Ln(X)

натуральный логарифм (X)

выражение

LnXP1(X)

натуральный логарифм(X+1)

Log10(X)

десятичный логарифм от X

Log2(X)

логарифм от Xпо основанию 2

LogN (N,X)

логарифм от Xпо основанию N

Max(A,B)

максимум двух чисел

int64, Single, Double

Extended

Min(A,B)

минимум двух чисел

Pi

число Пи

Poly(X,C)

вычисляет полином Xс массивом коэфф С

массив Double

Power (X, E)

Комментировать
589 просмотров
Комментариев нет, будьте первым кто его оставит

Это интересно
No Image Компьютеры
0 комментариев
No Image Компьютеры
0 комментариев
No Image Компьютеры
0 комментариев
No Image Компьютеры
0 комментариев