• Регистрация
hex oct bin
hex oct bin +17.80
н/д
  • Написать
  • Подписаться

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

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

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

Quake III Arena

    Теперь переместимся в уже далекий 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. Разберем это преобразование по подробнее. Для начала запишем число в экспоненциальной форме:

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

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

    Здесь L это разрядность мантиссы, а B смещение порядка, для F32 эти параметры будут 27-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. Вот в принципе и все - такая история.

Теги

    19.09.2020

    Комментарии

    • hex oct bin
      hex oct bin +17.80
      20.09.2020 11:37

      В редакторе все красиво: выравнивание, подсветка синтаксиса, жмешь опубликовать и все ломается

      • Sancho
        Sancho+105.07
        23.09.2020 17:31

        Передали разработчикам вопрос. Спасибо за интересный материал!

      • AlMich
        AlMich+24.00
        21.09.2020 13:41

        Здесь вот все лучше описано

        https://habr.com/ru/company/infopulse/blog/336110/

        еще комментарии надо прочитать