Алгоритм быстрого вычисления обратного квадратного корня и его актуальность

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

В ЭВМ вещественные числа представляются в формате IEEE754, согласно которому число представляется в виде бита знака, смещенного порядка и мантиссы. Для одинарной точности представления (соответствует типу float в языке C) справедливо следующее (рис. 1).

s	E – порядок	M – мантисса
1 бит	8 бит	23 бит

Рис. 1. Формат IEEE754

Формула для перехода из IEEE754 к десятичному формату:

M                       x-(-1)'-(1 + m )-2 e• где m = —, e = E - B, L = 223, B = 127 - смещение порядка.

Например, для десятичного числа 12.75 имеем:

12.7510 = 1100.112 = 1.10011 ∙ 23, тогда E = 3 + 127 = 13010 = 100000102, а M = 100112, т.е. получим (с учетом бита знака s = 0): 010000010100110000000000000000002.

Если интерпретировать порядок и мантиссу числа x как целые числа, то целочисленное представление числа x можно записать в виде

INT ( x ) = M + L · E .                              (2)

Пусть необходимо вычислить значение функции:

У = f ( x ) = x p , p = у , k ^e N .

Прологарифмируем обе части равенства по основанию 2:

log₂ y = log₂ x p = p • log₂ x .

Заменим x и y , используя выражение (1):

log 2 ( 1 + m_y ) 2 e y = p • log 2 ( 1 + m x ) 2 e x .

Упростим выражение, раскрыв логарифм произведения:

log 2 ( 1 ⁺ m y, ) ⁺ e y = p -[ l^og 2 ( 1 ^{+ m} _x ) ⁺ e x ] .                      ⁽³⁾

Рассмотрим функцию g ( m ) = log₂ ( 1 + m ) . В нашем случае мантисса нормализована, т.е. справедливо неравенство: 0 < m <  1 . На полуинтервале [0, 1) функция g ( m ) примерно ведет себя, как линейная функция m + ε. Меняя параметр ε, мы можем улучшить точность приближения, е = 0.0450465 дает лучшее приближение.

Таким образом, можно записать:

log 2 ( 1 ⁺ m x ) » m x ^+е .                             (4)

Преобразуем (3), используя (4):

e y ⁺ m y ^+е = p • ( e x ⁺ m x ^+е) .

Перейдем от e , m к E , M по равенствам, указанным в (1):

My            „ м,   )

E_v - B +—- + е = p - E_x - B + - - + е .

^y L ( x L J

Умножим обе части равенства на L и перенесем все слагаемые, не связанные с y , в правую часть равенства:

L - E y + M y = p - ( M x + L - E x ) + ( 1 - p ) - ( B -е ) - L .

Используя (2), получим:

INT ( y ) = p - INT ( x ) + ( 1 - p ) - ( B - e ) - L = p - INT ( x ) + const .      (5)

Значение |p| = -^ выбрано не просто так. Данный выбор позволяет заменить операцию умножения на сдвиг числа вправо или влево, в зависимости от знака p. Таким об- разом, целочисленную интерпретацию результата можно получить, сдвинув аргумент функции на k разрядов (вправо/влево) и прибавив константу. Данное приближение дает точность порядка 4%. Повысить точность можно с помощью метода Ньютона.

Рассчитаем значение константы из выражения (5) при p = —^- const =

²³ = 1597463007 10 = 0x5F3759DF .

Для уточнения результата используется метод Ньютона [2]. Формула для следующего приближения имеет следующий вид:

f ( Уп )

У п + 1 = У п - , x . ^{f (} У п )

Применим данный метод для функции f ( x ) = y = —^=;

x f (у, x )= y2—x=0;

f '( У , x ) =

—2

y 3 .

Подставив производную в (6), получим формулу для второго приближения результата:

2"— x

У п + 1 = У п ^— ^ п — Г = У п -( 1.5 ^— 0.5 x • У п² ) .

y n ³

Данное уточнение повышает точность результата до 0.18%.

Реализация алгоритма вычисления обратного квадратного корня на языке C:

inline float Q_Sqrt_Inverse(float x) { const float half = 0.5f * x; //запоминаем половину для последующей аппроксимации int y = * (int *) //приводим к int y = 0x5F3759DF - (y >> 1); //вычисляем первое приближение x = * (float *) //приводим к float x = x * (1.5f - (half * x * x)); //аппроксимация метода Ньютона return x;

}

В таблице 1 приведено полученное время выполнения функции по сравнению с функциями, использующими стандартные математические функции sqrt и sqrtf.

Время выполнения функций

Таблица 1

	sqrt	sqrtf	Q_Sqrt_Inverse
Без уточнения результата	0.170480	0.140246	0.060292
С уточнением результата	0.170480	0.140246	0.083915

В таблице 2 приведены коэффициенты ускорения приведенного алгоритма по отношению к стандартным функциям sqrt и sqrtf.

Таблица 2

	Ускорение относительно sqrt	Ускорение относительно sqrtf
Без уточнения результата	2.83	2.33
С уточнением результата	2.03	1.67

Коэффициенты ускорения алгоритма по отношению к стандартным функциям

Резюмируя полученные данные, можно сделать вывод, что в настоящее время, даже с учетом существования блоков FPU (Floating Point Unit), которые на аппаратном уровне реализуют большинство операций с плавающей точкой, приведенный алгоритм всё еще может использоваться по назначению. При точности 0.2% можно получить ускорение данной операции минимум на 50%.