Зведення числа в дійсну ступінь - delphi sources faq

Як, ніхто цього ще не придумав?
Є пропозиція
Чого ми досягли?
Апроксимація функції 2x
Новий варіант функції зведення в ступінь
Апроксимація функції log2x і "спеціалізація" піднесення до степеня
висновок
джерело мудрості

Як, ніхто цього ще не придумав?

Не беруся судити. Ймовірно, завдання про те, як максимально швидко звести дійсне додатне число в довільну дійсний ступінь, вирішувалася приблизно так само часто, як і вставала, - а вставала, гадаю, не раз. І все ж не так давно я з жахом виявив, що RTL зі складу Borland Delphi останніх версій (як Delphi 6, так і Delphi 7) підходить до вирішення не більше професійно, ніж старанний п'ятикласник, яка вперше зіштовхнулася з такою проблемою.

Погляньмо на вихідний код функції Power з модуля Math, люб'язно наданий Borland Software:

Примітно, що в благих цілях оптимізації процесор залишають наодинці з цілим натовпом розгалужень, що призводять його, врешті-решт, в загальному випадку до горезвісного рішенням п'ятикласника, а саме, до тривіальної формулою

Тут 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 * log2 (x + 1) для - (1-1 / sqrt (2))

Ось, загалом-то, і все. Але вже на перший погляд цього вистачає, щоб зрозуміти, що завдання може бути вирішена більш прямо, ніж пропонує RTL Borland Delphi.

Дійсно, чому не замінити показник ступеня в (1) на 2? Адже неперово число аж ніяк не є рідною для двійковій арифметики! тоді вийде

Цей вислів для x ** y відповідно до вищезазначеними п'ятьма інструкціями можна уявити як композицію функцій у такому вигляді:

Так як обчислити f (z) в одну дію неможливо, доводиться вважати так:

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

Перед виконанням цього фрагмента коду потрібно переконатися, що біти управління округленням в слові управління FPU встановлені в режим округлення до нуля. У Delphi це найпростіше зробити за допомогою функції SetRoundMode (модуль Math):

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

Повний код працездатною функції на Object Pascal виглядає так:

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

На жаль, виграш в швидкості абсолютно не відчувається. Це цілком зрозуміло: згідно з додатком 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) = an x n + an-1 x n-1 +. + A1 x + a0. Його коефіцієнти можуть бути постійні, а можуть деяким чином залежати від x. У першому випадку коефіцієнти легко знайти методом найменших квадратів, взявши значення вихідної функції в декількох точках і підібравши коефіцієнти так, щоб в цих точках многочлен брав значення, як можна ближчі до значень функції (докладний опис поліноміальною апроксимації функцій і методу найменших квадратів можна знайти в книгах, присвячених курсам обчислювальної математики або обробці експериментальних даних). Простота методу обертається істотним недоліком: він часом непоганий для виявлення якісних тенденцій, але погано відображає вихідну функцію кількісно, причому, як правило, похибка зростає зі зменшенням ступеня многочлена n. а швидкість обчислення g (x) з ростом n падає.

Гідна альтернатива, що дозволяє досить точно наближати гладкі криві, такі, як y = 2 ** x, - апроксимація сплайнами. Говорячи простою мовою (можливо, надто простим - нехай мене вибачать фахівці), сплайн - це крива, що моделює форму, прийняту пружним стержнем, деформованим шляхом закріплення в заданих точках. Вона проходить точно через задані точки, підкоряючись при цьому деяким додатковим умовам, зокрема, умові безперервності другої похідної. Існують різні види сплайнів. У цій роботі досить практично використання кубічних сплайнів. Кубічний сплайн на кожному відрізку між двома послідовними (в порядку зростання координати x) еталонними точками (x, f (x)) описується поліномом третього ступеня g (x) = a3 x 3 + a2 x 2 + a1 x + a0. де набір коефіцієнтів (a0, a1, a2, a3) свій для кожного відрізка. Пошук цих коефіцієнтів - не надто складне завдання, але опис методу її рішення виходить за рамки цієї статті. Таблиця коефіцієнтів. що виходить після обліку всіх зауважень цього розділу, додається до статті.

Отже, я обмежуся лише використанням отриманих мною значень коефіцієнтів. Щоб забезпечити необхідну точність на проміжку 0<=x <999, мне понадобились в качестве эталонных 2039 точек, которым соответствовали значения x=(i-40)/2. i= 0..2038. Сорок значений на отрицательной полуоси нужны были только для того, чтобы отразить поведение сплайна в этой части плоскости, слегка скорректировав таким образом его вид на остальных отрезках; в вычислениях эти 40 отрезков не участвуют, т.к. для значений x <0 можно воспользоваться (без ощутимого проигрыша в скорости или точности) соотношением 2**(-|x|)=1/(2**|x|).

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

Ось як це зробив я.

Як і для попередньої функції, потрібно забезпечити установку біт управління округленням в режим округлення до нуля.

Виконання цього фрагмента змінює вміст регістру EAX.

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

Малюнок 1. Максимальна похибка наближення функції Exp2 = 2 ** x (при x менше 990) не перевищує 0,004%.

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

Змінимо реалізацію зведення в ступінь відповідно до запропонованої аппроксимацией для 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).

Видно, що оцінена відносна похибка (фактично - відхилення від значення, що повертається вбудованої функцією) насправді не перевершує 0.004%!

У разі показника ступеня 17 коливання стають набагато частіше, однак загальна картина та сама.

Апроксимація функції log2x і "спеціалізація" піднесення до степеня

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

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

Якщо абсолютна величина x досить мала, члени ряду, вже починаючи з третього, досить слабо позначаються на результаті. Тому для значень x. досить близьких до 1 (щоб залишитися в обумовлених вище рамках прийнятних похибок, x повинен відстояти від 1 цієї статті не більше ніж на 0.01), обчислення log2 (x) = ln (x) / ln (2) = ln (x) * log2 (e ) = ln (1 + (x-1)) * log2 (e) можна замінити обчисленням (tt * t / 2) * log2 (e), де t = x-1.

Це дозволяє побудувати ще один варіант функції зведення в ступінь для значень підстави, близьких до 1. У ньому немає інструкції FYL2X, а замість неї присутній блок інструкцій, позначених символом "*" (знак "

"Означає наближене рівність):

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

Малюнок 4. Тимчасові витрати: New_Power_XNear1 - спеціалізований варіант New_Power.

На жаль, із зростанням показника ступеня максимальна похибка зростає, залишаючись, втім, в обумовлених межах (тобто менше 0,1%, більше того - менше 0,01%) навіть при дуже високих показниках:

Таким чином, нам вдалося отримати функції, переважаючі вбудовану за швидкістю від двох до чотирьох разів при похибки близько 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.