Ось, загалом-то, і все. Але вже на перший погляд цього
вистачає, щоб зрозуміти, що завдання може бути вирішена більш прямо, ніж пропонує
RTL Borland Delphi.
Дійсно, чому не замінити показник ступеня в
(1) на 2? Адже неперово число аж ніяк не є рідною для двійкової арифметики!
Тоді вийде
Це вираз для 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 виконується набагато швидше, ніж
на ранніх моделях, можна розраховувати, що навіть у тих ситуаціях, коли
доведеться перемежовувати виклик цього фрагмента коду з діями, які вимагають іншого
режиму округлення, при роботі на сучасній техніці це не призведе до
надмірним додатковим тимчасових витратах. Подробиці див., наприклад, в
[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.
Чого ми досягли?
Експерименти показали, що запропонований варіант
функції зведення до степеня підвищує точність обчислень на один-два знаки
після коми. Так як автору було кілька лінь писати повільний код для
надточної піднесення до степеня з метою перевірки точності запропонованого
алгоритму, то експеримент полягав у порівнянні результатів зі значенням,
що виходять в стандартному калькуляторі 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) = anxn + an-1xn-1 +...+ a1x + a0. Його
коефіцієнти можуть бути постійні, а можуть певним чином залежатиме від x. У
першому випадку коефіцієнти легко знайти методом найменших квадратів, взявши
значення вихідної функції в декількох точках і підібравши коефіцієнти так,
щоб у цих точках многочлен брав значення, як можна більш близькі до
значень функції (докладний опис поліноміальною апроксимації функцій і
методу найменших квадратів можна знайти в книгах, присвячених курсам
обчислювальної математики або обробці експериментальних даних). Простота
методу обертається суттєвий недолік: він часом непоганий для виявлення
якісних тенденцій, але погано відображає вихідну функцію кількісно,
причому, як правило, похибка зростає зі зменшенням ступеня многочлена n, а
швидкість обчислення g (x) з ростом n падає.
Гідна альтернатива, що дозволяє досить точно
наближати гладкі криві, такі, як y = 2 ** x, - апроксимація сплайнами. Говорячи
простою мовою (можливо, занадто простим - нехай мене вибачать фахівці),
сплайн - це крива, що моделює форму, що приймається пружним стержнем,
деформованим шляхом закріплення в заданих точках. Вона проходить точно через
задані точки, підкоряючись при цьому деяким додатковим умовам, в
Зокрема, умовою безперервності другої похідної. Існують різні види
сплайнів. У цій роботі досить практично використання кубічних сплайнів.
Кубічний сплайн на кожному відрізку між двома послідовними (у порядку
зростання координати 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 відрізків не беруть участі, тому що для значень x <0
можна скористатися (без відчутного програшу в швидкості чи точності)
співвідношенням 2 **(-| x |) = 1/(2 ** | x |).
Отже, у нас є таблиця коефіцієнтів, представлена
у вигляді масиву з 1999 блоків по 8 * 4 байт (якщо для подання коефіцієнтів
використовується тип 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 = abs (Exp2 (x) --
_Power (2, x))/_Power (2, x). Зрозуміло, вираз має сенс, якщо
_Power (2, x) <> 0.
Якщо побудувати графік відносної похибки,
стає видно, що в межах кожного з 1998 відрізків він має форму кривої
з одним максимумом, сходить до нуля на кінцях відрізка. При цьому межі
коливань величини похибки залишаються постійними на всіх відрізках, крім
декількох останніх - на них похибка зростає. Якщо не брати до
увагу ці відрізки, і обмежити область допустимих значень аргументу числом
990 (тобто x <990), то для опису поведінки відносної похибки в
залежно від x досить показати її графік на двох останніх допустимих для
значень 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;// abs (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) можна
замінити обчисленням (tt * t/2) * log2 (e), де t = x-1.
Це дозволяє побудувати ще один варіант функції
піднесення до степеня підстави для значень, близьких до 1. У ньому немає інструкції
FYL2X, а замість неї присутній блок інструкцій, позначених символом "*"
(знак "~" означає наближене рівність):
function New_Power_XNear1 (x, y: FLOATTYPE): FLOATTYPE;//
abs (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.
Висновок
Таким чином, нам вдалося отримати функції,
переважаючі вбудовану за швидкістю від двох до чотирьох разів при похибки
близько 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.
Для підготовки даної роботи були використані
матеріали з сайту http://www.rsdn.ru/
