Быстрый обратный корень

Сегодня мы поговорим про приближенные вычисления, Quake3 Arena, незаурядном инженере и программисте Джоне Кармаке. И даже вспомним одного из разработчиков MATLAB.

Итак, тема нашего нашего разговора - быстрый способ приближенного вычисления X^-0.5. Во многих прикладных областях нормировка на квадратный корень из какой-либо величины используется довольно часто. Например, про обработки сигналов и говорить нечего, возьмем корреляцию, которая встречается практически повсеместно. На что ее нормировать? На квадратный корень из произведения квадратичного произведения сигналов (Айфечер, 5я глава, формула 5.8). В компьютерной графике дела обстоят схожим образом, упомянем на всякий случай поверхность нормалей.

Теперь переместимся в уже далекий 1999 год, а точнее самое начало 2000х. Времена компьютерных клубов, "кваки" и "каунтер стрике 1.1". Когда можно было без проблем купить "отвертку" для папы, а главной целью в жизни было отпроситься у мамки на ночь в компьютерный клуб. Скажу честно, мне это так и не удалось, сходил на ночь я уже на свои и та магическая аура, исходившая от компьютеров Pentium III, уже улетучилась. 

 В 1999 году Id Software выпускают прорывной шутер, давший невероятный толчок развитию не только игр, но и компьютерной графики в целом. В реализации движка использовался тот самый быстрый обратный корень. Когда реализация игры была выложена в открытый доступ, то около данной функции был обнаружен комментарий, содержащий обсценное слово на английском языке, обозначающее половое сношение. Очевидно, Джон Кармак не разобрался в реализации, но она работала и работала хорошо. Про Кармака советую почитать отдельно, уж очень он незаурядная личность.

А теперь сам алгоритм. Это тот самый код из Quake III Арена и те самые комментарии, во многих русскоязычных источниках их удаляют, будем надеяться, у администрации платформы тестикулы на месте.

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

В принципе на этом статью можно заканчивать, алгоритм перед вами и он работает, сам проверял. Но все-таки давайте разберемся, что же тут происходит. 

Во-первых, это отображение числа с плавающей точкой, как 32х битное целое число. В оригинальном алгоритме это сделано довольно витьевато, я предпочел бы применить union.

union {
		float f;
		int i;
	} conv = {number};

Признаюсь честно, для меня стало открытием, что такое преобразование дает довольно неплохую аппроксимацию логарифма числа по основанию 2. Разберем это преобразование по подробнее. Для начала запишем число в экспоненциальной форме:

Здесь m_x мантисса а e_x соответственно экспонента. Отметим также, что знаковый бит всегда должен быть равен нулю. Тогда логарифм от такого числа можно представить, как

Тогда биекцию между целочисленным представлением и представлением с плавающей точкой можно представить в следующем виде:

Здесь L это разрядность мантиссы, а B смещение порядка, для F32 эти параметры будут 2⁷-1 соответственно. Также мы уже два раза встречается сигма, это не что иное, как параметр приближения. Теперь выразим логарифм по основанию 2 от x.

Все выражения, записанные выше приведены для y=x, давайте теперь подставим в них

y=x^-0.5.

Этой формулой, как раз и описывается та самая строчка:

i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 

Точное происхождение именно это константы, установить так и не удалось. Далее проведем обратное преобразование:

В коде Кармака это выглядит следующим образом:

y  = * ( float * ) &i;

Затем проводится одна итерация приближения по методу Ньютона:

y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

Метод широко распространен и описан во множестве источников, возможно, мы разберем его в будущем.

На рисунке выше приводится ошибка в расчетах с использованием приведенного метода. 

Стоит только отметить, что многие ошибочно приписывают разработку этого алгоритма самому Кармаку, однако это неверно. Сам Джон предполагал, что алгоритм принес в компанию Майкл Абраш. Самая ранняя обнаруженная версия реализации алгоритма написана Гэрри Таролли, но также есть предположение, что алгоритм разработали его коллеги по Ardent Computer Грег Уолш и Клив Моулер. Кстати, Клив Моулер в последствии разработает тот самый MATLAB. Вот в принципе и все - такая история.