В этой статье мы поговорим о «магической» константе 0x5f3759df, лежащей в основе элегантного алгоритмического трюка для быстрого вычисления обратного квадратного корня.
Вот полная реализация этого алгоритма:
Этот код вычисляет некоторое (достаточно неплохое) приближение для формулы
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/460/bc7/d55/460bc7d55c4920d2de99f2086fdd06b9.png)
Сегодня данная реализация уже хорошо известна, и стала она такой после появления в коде игры Quake III Arena в 2005 году. Её создание когда-то приписывали Джону Кармаку, но выяснилось, что корни уходят намного дальше – к Ardent Computer, где в середине 80-ых её написал Грег Уолш. Конкретно та версия кода, которая показана выше (с забавными комментариями), действительно из кода Quake.
В этой статье мы попробуем разобраться с данным хаком, математически вывести эту самую константу и попробовать обобщить данный метод для вычисления произвольных степеней от -1 до 1.
Да, понадобится немного математики, но школьного курса будет более, чем достаточно.
Зачем вообще нам может понадобиться считать обратный квадратный корень, да ещё и пытаться делать это настолько быстро, что нужно реализовывать для этого специальные хаки? Потому, что это одна из основных операций в 3D программировании. При работе с 3D графикой используются нормали поверхностей. Это вектор (с тремя координатами) длиной в единицу, который нужен для описания освещения, отражения и т.д. Таких векторов нужно много. Нет, даже не просто много, а МНОГО. Как мы нормализуем (приводим длину к единице) вектор? Мы делим каждую координату на текущую длину вектора. Ну или, перефразируя, умножаем каждую координату вектора на величину:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/aba/d41/552/abad41552fcd31d1b1886b047aafc24a.png)
Расчёт
относительно прост и работает быстро на всех типах процессоров. А вот расчёт квадратного корня и деление – дорогие операции. И вот поэтому – встречайте алгоритм быстрого извлечения обратного квадратного корня — FastInvSqrt.
Что же делает вышеуказанная функция для вычисления результата? Она состоит из четырёх основных шагов. Первым делом она берёт входное число (которое пришло нам в формате float) и интерпретирует его биты как значение новой переменной i типа integer (целое число).
Далее над полученным целым числом выполняется некоторая целочисленная арифметическая операция, которая работает достаточно быстро и даёт нам некоторую аппроксимацию требуемого результата
То, что мы получили в результате данной операции, ещё не является, собственно, результатом. Это лишь целое число, биты которого представляют некоторое другое число с плавающей запятой, которое нам нужно. А значит, необходимо выполнить обратное преобразование int во float.
И, наконец, выполняется одна итерация метода Ньютона для улучшения аппроксимации.
И вот теперь у нас имеется отличная аппроксимация операции извлечения обратного квадратного корня. Последняя часть алгоритма (метод Ньютона) – достаточно тривиальная вещь, я не буду на этом останавливаться. Ключевой частью алгоритма является шаг №2: выполнение некоторой хитрой арифметической операции над целым числом, полученным от интерпретации битов числа с плавающей запятой в качестве содержимого типа int. На этом мы и сфокусируемся.
Для понимания дальнейшего текста нам нужно вспомнить формат, в котором хранятся в памяти числа с плавающей запятой. Я опишу здесь только то, что важно здесь для нас (остальное всегда можно подсмотреть в Википедии). Число с плавающей запятой хранится как комбинация трёх составляющих: знака, экспоненты и мантиссы. Вот биты 32-битного числа с плавающей запятой:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/749/cc6/41e/749cc641eb4d5dafd085e8c23f8826aa.png)
Первым идёт бит знака, дальше 8 битов экспоненты и 23 бита мантиссы. Поскольку мы здесь имеем дело с вычислением квадратного корня, то я предположу, что мы будем иметь дело лишь с положительными числами (первый бит всегда будет 0).
Рассматривая число с плавающей запятой как просто набор битов, экспонента и мантисса могут восприниматься как просто два положительных целых числа. Давайте обозначим их, соответственно, Е и М (поскольку дальше мы будем часто на них ссылаться). С другой стороны, интерпретируя биты числа с плавающей запятой, мы будем рассматривать мантиссу как число между 0 и 1, т.е. все нули в мантиссе будут означать 0, а все единицы – некоторое число очень близкое (но всё же не равное) 1. Ну и вместо того, чтобы интерпретировать экспоненту как беззнаковое 8-битное целое число, давайте вычтем смещение (обозначим его как В), чтобы получить знаковое целое число в диапазоне от -127 до 128. Давайте обозначим float-интерпретацию этих значений как е и m. Чтобы не путаться, будем использовать прописные обозначения (Е, М) для обозначения целочисленных значений и строчные (е, m) – для обозначения чисел с плавающей запятой (float).
Преобразование из одного в другое тривиально:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/e20/4f6/0da/e204f60da8cd0817a4e780d1bb040767.png)
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/b3c/587/908/b3c5879085ad90e7bdb51854803b2375.png)
В этих формулах для 32-битных чисел L = 2^23, а В = 127. Имея некоторые значения e и m, можно получить само представляемое ими число:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/1e4/d11/3b3/1e4d113b311c654d1aa11bf6ab69aad9.png)
и значение соответствующей им целочисленной интерпретации числа:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/e7f/4fb/7b4/e7f4fb7b4e0d26ea4fa5f0c87484a881.png)
Теперь у нас есть почти все кусочки паззла, которые нужны для объяснения “хака” в коде выше. Итак, нам на вход приходит некоторое число х и требуется рассчитать его обратный квадратный корень:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/fc5/217/c7e/fc5217c7ea90261a96111ab15890d1e5.png)
По некоторым причинам, которые вскоре станут понятны, я начну с того, что возьму логарифм по основанию 2 от обеих частей данного уравнения:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/2e8/142/d4a/2e8142d4a0691814fa7e1fb989e2c0f9.png)
Поскольку числа, с которыми мы работаем, на самом деле являются числами с плавающей запятой, мы можем представить х и у согласно вышеуказанной формуле представления таких чисел:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/e3e/483/0a1/e3e4830a1140ec9e98526d80b71b7fac.png)
Ох уж эти логарифмы. Возможно, вы не пользовались ими со школьных времён и немного подзабыли. Не волнуйтесь, немного дальше мы избавимся от них, но пока что они нам ещё нужны, так что давайте посмотрим, как они здесь работают.
В обеих частях уравнения у нас находится выражение вида:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/7f0/71a/3cb/7f071a3cb749faaa906b062c6573685b.png)
где v находится в диапазоне от 0 до 1. Можно заметить, что для v от 0 до 1 эта функция достаточно близка к прямой линии:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/230/ec1/a1e/230ec1a1e4f2af4784e816add3307403.png)
Ну или в виде выражения:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/d62/ce4/ee7/d62ce4ee77931f98c5d6ab2f6c29a68e.png)
Где σ – некоторая константа. Это не идеальное приближение, но мы можем попытаться подобрать σ таким образом, чтобы оно было достаточно неплохим. Теперь, использовав его, мы можем преобразовать вышеуказанное равенство с логарифмами в другое, уже не строго равное, но достаточно близкое и, самое главное, СТРОГО ЛИНЕЙНОЕ выражение:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/81e/439/481/81e439481fa38dd0eee19213002c8081.png)
Это уже что-то! Теперь самое время прекратить работать с представлениями в виде чисел с плавающей запятой и перейти к целочисленному представлению мантиссы и экспоненты:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/802/de1/219/802de1219ff7f298350ffb67675b3051.png)
Выполнив ещё несколько тривиальных преобразований (можно пропустить детали), мы получим нечто, выглядящее уже довольно знакомо:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/ae8/7c6/387/ae87c6387c6e08b48666c9f2ca86665a.png)
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/52c/844/4ed/52c8444ed66d49d544c7ec05734de8e2.png)
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/7fa/f54/acf/7faf54acfd2548eeb99388992ab54f32.png)
Посмотрите внимательно на левую и правую части последнего уравнения. Как мы видим, у нас получилось выражение целочисленной формы требуемого значения y, выраженное через линейного вида формулу, включающую целочисленное представление значения х:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/9fc/897/ff2/9fc897ff2e6b16793c55b3e3df569c9c.png)
Говоря простыми словами: “y (в целочисленной форме) – это некоторая константа минус половина от целочисленной формы х”. В виде кода это:
Очень похоже на формулу в коде функции в начале статьи, правда?
Нам осталось найти константу K. Мы уже знаем значения B и L, но ещё не знаем чему равно σ.
Как вы помните, σ – это некоторое “поправочное значение”, которое мы ввели для улучшения аппроксимации функции логарифма к прямой линии на отрезке от 0 до 1. Т.е. мы можем подобрать это число сами. Я возьму число 0.0450465, как дающее неплохое приближение и использованное в оригинальной реализации. Используя его, мы получим:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/245/b6d/973/245b6d973e775793897d08295d28517e.png)
Угадаете, как число 1597463007 представляется в HEX? Ну, конечно, это 0x5f3759df. Ну, так и должно было получиться, поскольку я выбрал σ таким образом, чтобы получить именно это число.
Таким образом, данное число не является битовой маской (как думают некоторые люди просто потому, что оно записано в hex-форме), а результатом вычисления аппроксимации.
Но, как сказал бы Кнут: “Мы пока что лишь доказали, что это должно работать, но не проверили, что это действительно работает”. Чтобы оценить качество нашей формулы, давайте нарисуем графики вычисленного таким образом значения обратного квадратного корня и настоящей, точной его реализации:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/92f/34d/b69/92f34db69cf09309e4bb98a243babe43.png)
График построен для чисел от 1 до 100. Неплохо, правда? И это никакая не магия, не фокус, а просто правильное использование несколько, возможно, экзотичного трюка с представлением чисел с плавающей запятой в виде целочисленных и наоборот.
Все вышеуказанные преобразования и выражения дали нам не только объяснение константы 0x5f3759df, но и ещё несколько ценных выводов.
Во-первых, поговорим о значениях чисел L и B. Они определяются не нашей задачей по извлечению обратного квадратного корня, а форматом хранения числа с плавающей запятой. Это означает, что тот же трюк может быть проделан и для 64-битных и для 128-битных чисел с плавающей запятой – нужно лишь повторить вычисления для рассчета других констант.
Во-вторых, нам не очень-то и важно выбранное значение σ. Оно может не давать (да и на самом деле – не даёт) лучшую аппроксимацию функции x + σ к логарифму. σ выбрано таким, поскольку оно даёт лучший результат совместно с последующим применением алгоритма Ньютона. Если бы мы его не применяли, то выбор σ являлся бы отдельной интересной задачей сам по себе, эта тема раскрыта в других публикациях.
Ну и в конце-концов, давайте посмотрим на коэффициент “-1/2” в финальной формуле. Он получился таким из-за сути того, что мы хотели вычислить (“обратного квадратного корня”). Но, в общем, степень здесь может быть любой от -1 до 1. Если мы обозначим степень как р и обобщим все те же самые преобразования, то вместо “-1/2” мы получим:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/60b/dde/096/60bdde0967db98dcbbef0cf0ae9c73c8.png)
Давайте положим р=0.5 Это будет вычисление обычного (не обратного) квадратного корня числа:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/327/8ac/3fe/3278ac3fe4117548e0d07bbd7a3cb043.png)
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/1b0/722/def/1b0722def05d4e8ed10e8c2f34228e13.png)
В виде кода это будет:
И что, это работает? Конечно, работает:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/161/16b/aec/16116baec153f85a3445a22fe1fa2cc3.png)
Это, наверное, хорошо известный способ быстрой аппроксимации значения квадратного корня, но беглый гуглинг не дал мне его названия. Возможно, вы подскажете?
Данный способ будет работать и с более “странными” степенями, вроде кубического корня:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/3cd/7e0/5c3/3cd7e05c3d16a8f90d7f75c99e5888df.png)
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/ddf/714/c2a/ddf714c2a23d6abe0ca4b77e0ae5c7d0.png)
Что в коде будет выражено, как:
К сожалению, из-за степени 1/3 мы не можем использовать битовые операции сдвига и вынуждены применить здесь умножение на 0.333f. Аппроксимация всё ещё достаточно хороша:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/237/28c/c8f/23728cc8f0790f93860291de889d6a50.png)
К этому моменту вы уже могли заменить, что изменение степени вычисляемой функции достаточно тривиально меняет наши вычисления: мы просто рассчитываем новую константу. Это совершенно не затратная операция и мы вполне можем делать это даже на этапе выполнения кода, для разных требуемых степеней. Если мы перемножим две заранее известных константы:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/b60/418/678/b60418678d4b84e48afdaf07d852cf3a.png)
То сможем вычислять требуемые значения на лету, для произвольной степени от -1 до 1:
Чуть упростив выражение, мы даже можем сэкономить на одном умножении:
Эта формула даёт нам “магическую” константу, с помощью которой можно рассчитывать на лету разные степени чисел (для степеней от -1 до 1). Для полного счастья нам не хватает лишь уверенности в том, что вычисленное таким образом приближенное значение может быть столь же эффективно улучшено алгоритмом Ньютона, как это происходило в оригинальной реализации для обратного квадратного корня. Я не изучал эту тему более глубоко и это, вероятно, будет темой отдельной публикации (наверное, не моей).
Выражение выше содержит новую “магическую константу” 0x3f7a3bea. В некотором смысле (из-за своей универсальности) она даже “более магическая”, чем константа в оригинальном коде. Давайте назовём её С и посмотрим на неё чуть внимательнее.
Давайте проверим работу нашей формулы для случая, когда p = 0. Как вы помните из курса математики, любое число в степени 0 равно единице. Что же будет с нашей формулой? Всё очень просто – умножение на 0 уничтожит второе слагаемое и у нас останется:
Что, действительно является константой и, будучи переведённой во float-формат, даст нам 0.977477 – т.е. “почти 1”. Поскольку мы имеем дело с аппроксимациями – это неплохое приближение. Кроме того, это говорит нам кое-что ещё. Наша константа С имеет совершенно не случайное значение. Это единица в формате чисел с плавающей запятой (ну или “почти единица”)
Это интересно. Давайте взглянем поближе:
Целочисленное представление C это:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/8f8/04f/625/8f804f625d5a14378881ce21f12d1129.png)
Это почти, но всё же не совсем, форма числа с плавающей запятой. Единственная проблема – мы вычитаем вторую часть выражения, а должны были бы её прибавлять. Но это можно исправить:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/54d/dc1/4e4/54ddc14e43de0edefb16a7e8fcef704f.png)
Вот теперь это выглядит в точности как целочисленное представление числа с плавающей запятой. Чтобы определить, какого именно числа, мы вычислим экспоненту и мантиссу, а потом уже и само число С. Вот экспонента:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/181/499/46d/18149946d0b62fd204ecf880a65c6744.png)
А вот мантисса:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/75f/f7d/c02/75ff7dc02e9cea442e4fa32953a9852b.png)
А, значит, значение самого числа будет равно:
![image](https://habrastorage.org/r/w1560/getpro/habr/post_images/1cb/c9a/34a/1cbc9a34acef9e94c7e349ffb35ece89.png)
И правда, если поделить наше σ (а оно было равно 0.0450465) на 2 и отнять результат от единицы, то мы получим уже известное нам 0.97747675, то самое, которое “почти 1”. Это позволяет нам посмотреть на С с другой стороны и вычислять его на рантайме:
Учтите, что для фиксированного σ все эти числа будут константами и компилятор сможет рассчитать их на этапе компиляции. Результатом будет 0x3f7a3beb, что не совсем 0x3f7a3bea из расчетов выше, но отличается от него всего на 1 бит (наименее значимый). Получить из этого числа оригинальную константу из кода (и заголовка данной статьи) можно, умножив результат вычислений ещё на 1.5.
Со всеми этими выкладками мы приблизились к пониманию того, что в коде в начале статьи нет никакой “магии”, “фокусов”, “интуиции”, “подбора” и прочих грязных хаков, а есть лишь чистая математика, во всей своей первозданной красоте. Для меня главным выводом из этой истории стала новость о том, что преобразование float в int и наоборот путем переинтерпретации одного и того же набора бит – это не ошибка программиста и не “хак”, а вполне себе иногда разумная операция. Слегка, конечно, экзотическая, но зато очень быстрая и дающая практически полезные результаты. И, мне кажется, для этой операции могут быть найдены и другие применения – будем ждать.
Вот полная реализация этого алгоритма:
float FastInvSqrt(float x) {
float xhalf = 0.5f * x;
int i = *(int*)&x; // представим биты float в виде целого числа
i = 0x5f3759df - (i >> 1); // какого черта здесь происходит ?
x = *(float*)&i;
x = x*(1.5f-(xhalf*x*x));
return x;
}
Этот код вычисляет некоторое (достаточно неплохое) приближение для формулы
![image](https://habrastorage.org/getpro/habr/post_images/460/bc7/d55/460bc7d55c4920d2de99f2086fdd06b9.png)
Сегодня данная реализация уже хорошо известна, и стала она такой после появления в коде игры Quake III Arena в 2005 году. Её создание когда-то приписывали Джону Кармаку, но выяснилось, что корни уходят намного дальше – к Ardent Computer, где в середине 80-ых её написал Грег Уолш. Конкретно та версия кода, которая показана выше (с забавными комментариями), действительно из кода Quake.
В этой статье мы попробуем разобраться с данным хаком, математически вывести эту самую константу и попробовать обобщить данный метод для вычисления произвольных степеней от -1 до 1.
Да, понадобится немного математики, но школьного курса будет более, чем достаточно.
Зачем?
Зачем вообще нам может понадобиться считать обратный квадратный корень, да ещё и пытаться делать это настолько быстро, что нужно реализовывать для этого специальные хаки? Потому, что это одна из основных операций в 3D программировании. При работе с 3D графикой используются нормали поверхностей. Это вектор (с тремя координатами) длиной в единицу, который нужен для описания освещения, отражения и т.д. Таких векторов нужно много. Нет, даже не просто много, а МНОГО. Как мы нормализуем (приводим длину к единице) вектор? Мы делим каждую координату на текущую длину вектора. Ну или, перефразируя, умножаем каждую координату вектора на величину:
![image](https://habrastorage.org/getpro/habr/post_images/aba/d41/552/abad41552fcd31d1b1886b047aafc24a.png)
Расчёт
![image](https://habrastorage.org/getpro/habr/post_images/ac0/5bb/8b4/ac05bb8b432b392343d3cb6797c929ce.png)
Что он делает?
Что же делает вышеуказанная функция для вычисления результата? Она состоит из четырёх основных шагов. Первым делом она берёт входное число (которое пришло нам в формате float) и интерпретирует его биты как значение новой переменной i типа integer (целое число).
int i = *(int*)&x; // представим биты float в виде целого числа
Далее над полученным целым числом выполняется некоторая целочисленная арифметическая операция, которая работает достаточно быстро и даёт нам некоторую аппроксимацию требуемого результата
i = 0x5f3759df - (i >> 1); // какого черта здесь происходит ?
То, что мы получили в результате данной операции, ещё не является, собственно, результатом. Это лишь целое число, биты которого представляют некоторое другое число с плавающей запятой, которое нам нужно. А значит, необходимо выполнить обратное преобразование int во float.
x = *(float*)&i;
И, наконец, выполняется одна итерация метода Ньютона для улучшения аппроксимации.
x = x*(1.5f-(xhalf*x*x));
И вот теперь у нас имеется отличная аппроксимация операции извлечения обратного квадратного корня. Последняя часть алгоритма (метод Ньютона) – достаточно тривиальная вещь, я не буду на этом останавливаться. Ключевой частью алгоритма является шаг №2: выполнение некоторой хитрой арифметической операции над целым числом, полученным от интерпретации битов числа с плавающей запятой в качестве содержимого типа int. На этом мы и сфокусируемся.
Какого черта здесь происходит?
Для понимания дальнейшего текста нам нужно вспомнить формат, в котором хранятся в памяти числа с плавающей запятой. Я опишу здесь только то, что важно здесь для нас (остальное всегда можно подсмотреть в Википедии). Число с плавающей запятой хранится как комбинация трёх составляющих: знака, экспоненты и мантиссы. Вот биты 32-битного числа с плавающей запятой:
![image](https://habrastorage.org/getpro/habr/post_images/749/cc6/41e/749cc641eb4d5dafd085e8c23f8826aa.png)
Первым идёт бит знака, дальше 8 битов экспоненты и 23 бита мантиссы. Поскольку мы здесь имеем дело с вычислением квадратного корня, то я предположу, что мы будем иметь дело лишь с положительными числами (первый бит всегда будет 0).
Рассматривая число с плавающей запятой как просто набор битов, экспонента и мантисса могут восприниматься как просто два положительных целых числа. Давайте обозначим их, соответственно, Е и М (поскольку дальше мы будем часто на них ссылаться). С другой стороны, интерпретируя биты числа с плавающей запятой, мы будем рассматривать мантиссу как число между 0 и 1, т.е. все нули в мантиссе будут означать 0, а все единицы – некоторое число очень близкое (но всё же не равное) 1. Ну и вместо того, чтобы интерпретировать экспоненту как беззнаковое 8-битное целое число, давайте вычтем смещение (обозначим его как В), чтобы получить знаковое целое число в диапазоне от -127 до 128. Давайте обозначим float-интерпретацию этих значений как е и m. Чтобы не путаться, будем использовать прописные обозначения (Е, М) для обозначения целочисленных значений и строчные (е, m) – для обозначения чисел с плавающей запятой (float).
Преобразование из одного в другое тривиально:
![image](https://habrastorage.org/getpro/habr/post_images/e20/4f6/0da/e204f60da8cd0817a4e780d1bb040767.png)
![image](https://habrastorage.org/getpro/habr/post_images/b3c/587/908/b3c5879085ad90e7bdb51854803b2375.png)
В этих формулах для 32-битных чисел L = 2^23, а В = 127. Имея некоторые значения e и m, можно получить само представляемое ими число:
![image](https://habrastorage.org/getpro/habr/post_images/1e4/d11/3b3/1e4d113b311c654d1aa11bf6ab69aad9.png)
и значение соответствующей им целочисленной интерпретации числа:
![image](https://habrastorage.org/getpro/habr/post_images/e7f/4fb/7b4/e7f4fb7b4e0d26ea4fa5f0c87484a881.png)
Теперь у нас есть почти все кусочки паззла, которые нужны для объяснения “хака” в коде выше. Итак, нам на вход приходит некоторое число х и требуется рассчитать его обратный квадратный корень:
![image](https://habrastorage.org/getpro/habr/post_images/fc5/217/c7e/fc5217c7ea90261a96111ab15890d1e5.png)
По некоторым причинам, которые вскоре станут понятны, я начну с того, что возьму логарифм по основанию 2 от обеих частей данного уравнения:
![image](https://habrastorage.org/getpro/habr/post_images/2e8/142/d4a/2e8142d4a0691814fa7e1fb989e2c0f9.png)
Поскольку числа, с которыми мы работаем, на самом деле являются числами с плавающей запятой, мы можем представить х и у согласно вышеуказанной формуле представления таких чисел:
![image](https://habrastorage.org/getpro/habr/post_images/e3e/483/0a1/e3e4830a1140ec9e98526d80b71b7fac.png)
Ох уж эти логарифмы. Возможно, вы не пользовались ими со школьных времён и немного подзабыли. Не волнуйтесь, немного дальше мы избавимся от них, но пока что они нам ещё нужны, так что давайте посмотрим, как они здесь работают.
В обеих частях уравнения у нас находится выражение вида:
![image](https://habrastorage.org/getpro/habr/post_images/7f0/71a/3cb/7f071a3cb749faaa906b062c6573685b.png)
где v находится в диапазоне от 0 до 1. Можно заметить, что для v от 0 до 1 эта функция достаточно близка к прямой линии:
![image](https://habrastorage.org/getpro/habr/post_images/230/ec1/a1e/230ec1a1e4f2af4784e816add3307403.png)
Ну или в виде выражения:
![image](https://habrastorage.org/getpro/habr/post_images/d62/ce4/ee7/d62ce4ee77931f98c5d6ab2f6c29a68e.png)
Где σ – некоторая константа. Это не идеальное приближение, но мы можем попытаться подобрать σ таким образом, чтобы оно было достаточно неплохим. Теперь, использовав его, мы можем преобразовать вышеуказанное равенство с логарифмами в другое, уже не строго равное, но достаточно близкое и, самое главное, СТРОГО ЛИНЕЙНОЕ выражение:
![image](https://habrastorage.org/getpro/habr/post_images/81e/439/481/81e439481fa38dd0eee19213002c8081.png)
Это уже что-то! Теперь самое время прекратить работать с представлениями в виде чисел с плавающей запятой и перейти к целочисленному представлению мантиссы и экспоненты:
![image](https://habrastorage.org/getpro/habr/post_images/802/de1/219/802de1219ff7f298350ffb67675b3051.png)
Выполнив ещё несколько тривиальных преобразований (можно пропустить детали), мы получим нечто, выглядящее уже довольно знакомо:
![image](https://habrastorage.org/getpro/habr/post_images/ae8/7c6/387/ae87c6387c6e08b48666c9f2ca86665a.png)
![image](https://habrastorage.org/getpro/habr/post_images/52c/844/4ed/52c8444ed66d49d544c7ec05734de8e2.png)
![image](https://habrastorage.org/getpro/habr/post_images/7fa/f54/acf/7faf54acfd2548eeb99388992ab54f32.png)
Посмотрите внимательно на левую и правую части последнего уравнения. Как мы видим, у нас получилось выражение целочисленной формы требуемого значения y, выраженное через линейного вида формулу, включающую целочисленное представление значения х:
![image](https://habrastorage.org/getpro/habr/post_images/9fc/897/ff2/9fc897ff2e6b16793c55b3e3df569c9c.png)
Говоря простыми словами: “y (в целочисленной форме) – это некоторая константа минус половина от целочисленной формы х”. В виде кода это:
i = K - (i >> 1);
Очень похоже на формулу в коде функции в начале статьи, правда?
Нам осталось найти константу K. Мы уже знаем значения B и L, но ещё не знаем чему равно σ.
Как вы помните, σ – это некоторое “поправочное значение”, которое мы ввели для улучшения аппроксимации функции логарифма к прямой линии на отрезке от 0 до 1. Т.е. мы можем подобрать это число сами. Я возьму число 0.0450465, как дающее неплохое приближение и использованное в оригинальной реализации. Используя его, мы получим:
![image](https://habrastorage.org/getpro/habr/post_images/245/b6d/973/245b6d973e775793897d08295d28517e.png)
Угадаете, как число 1597463007 представляется в HEX? Ну, конечно, это 0x5f3759df. Ну, так и должно было получиться, поскольку я выбрал σ таким образом, чтобы получить именно это число.
Таким образом, данное число не является битовой маской (как думают некоторые люди просто потому, что оно записано в hex-форме), а результатом вычисления аппроксимации.
Но, как сказал бы Кнут: “Мы пока что лишь доказали, что это должно работать, но не проверили, что это действительно работает”. Чтобы оценить качество нашей формулы, давайте нарисуем графики вычисленного таким образом значения обратного квадратного корня и настоящей, точной его реализации:
![image](https://habrastorage.org/getpro/habr/post_images/92f/34d/b69/92f34db69cf09309e4bb98a243babe43.png)
График построен для чисел от 1 до 100. Неплохо, правда? И это никакая не магия, не фокус, а просто правильное использование несколько, возможно, экзотичного трюка с представлением чисел с плавающей запятой в виде целочисленных и наоборот.
Но и это ещё не всё!
Все вышеуказанные преобразования и выражения дали нам не только объяснение константы 0x5f3759df, но и ещё несколько ценных выводов.
Во-первых, поговорим о значениях чисел L и B. Они определяются не нашей задачей по извлечению обратного квадратного корня, а форматом хранения числа с плавающей запятой. Это означает, что тот же трюк может быть проделан и для 64-битных и для 128-битных чисел с плавающей запятой – нужно лишь повторить вычисления для рассчета других констант.
Во-вторых, нам не очень-то и важно выбранное значение σ. Оно может не давать (да и на самом деле – не даёт) лучшую аппроксимацию функции x + σ к логарифму. σ выбрано таким, поскольку оно даёт лучший результат совместно с последующим применением алгоритма Ньютона. Если бы мы его не применяли, то выбор σ являлся бы отдельной интересной задачей сам по себе, эта тема раскрыта в других публикациях.
Ну и в конце-концов, давайте посмотрим на коэффициент “-1/2” в финальной формуле. Он получился таким из-за сути того, что мы хотели вычислить (“обратного квадратного корня”). Но, в общем, степень здесь может быть любой от -1 до 1. Если мы обозначим степень как р и обобщим все те же самые преобразования, то вместо “-1/2” мы получим:
![image](https://habrastorage.org/getpro/habr/post_images/60b/dde/096/60bdde0967db98dcbbef0cf0ae9c73c8.png)
Давайте положим р=0.5 Это будет вычисление обычного (не обратного) квадратного корня числа:
![image](https://habrastorage.org/getpro/habr/post_images/327/8ac/3fe/3278ac3fe4117548e0d07bbd7a3cb043.png)
![image](https://habrastorage.org/getpro/habr/post_images/1b0/722/def/1b0722def05d4e8ed10e8c2f34228e13.png)
В виде кода это будет:
i = 0x1fbd1df5 + (i >> 1);
И что, это работает? Конечно, работает:
![image](https://habrastorage.org/getpro/habr/post_images/161/16b/aec/16116baec153f85a3445a22fe1fa2cc3.png)
Это, наверное, хорошо известный способ быстрой аппроксимации значения квадратного корня, но беглый гуглинг не дал мне его названия. Возможно, вы подскажете?
Данный способ будет работать и с более “странными” степенями, вроде кубического корня:
![image](https://habrastorage.org/getpro/habr/post_images/3cd/7e0/5c3/3cd7e05c3d16a8f90d7f75c99e5888df.png)
![image](https://habrastorage.org/getpro/habr/post_images/ddf/714/c2a/ddf714c2a23d6abe0ca4b77e0ae5c7d0.png)
Что в коде будет выражено, как:
i = (int) (0x2a517d3c + (0.333f * i));
К сожалению, из-за степени 1/3 мы не можем использовать битовые операции сдвига и вынуждены применить здесь умножение на 0.333f. Аппроксимация всё ещё достаточно хороша:
![image](https://habrastorage.org/getpro/habr/post_images/237/28c/c8f/23728cc8f0790f93860291de889d6a50.png)
И даже более того!
К этому моменту вы уже могли заменить, что изменение степени вычисляемой функции достаточно тривиально меняет наши вычисления: мы просто рассчитываем новую константу. Это совершенно не затратная операция и мы вполне можем делать это даже на этапе выполнения кода, для разных требуемых степеней. Если мы перемножим две заранее известных константы:
![image](https://habrastorage.org/getpro/habr/post_images/b60/418/678/b60418678d4b84e48afdaf07d852cf3a.png)
То сможем вычислять требуемые значения на лету, для произвольной степени от -1 до 1:
i = (1 - p) * 0x3f7a3bea + (p * i);
Чуть упростив выражение, мы даже можем сэкономить на одном умножении:
i = 0x3f7a3bea + p * (i - 0x3f7a3bea);
Эта формула даёт нам “магическую” константу, с помощью которой можно рассчитывать на лету разные степени чисел (для степеней от -1 до 1). Для полного счастья нам не хватает лишь уверенности в том, что вычисленное таким образом приближенное значение может быть столь же эффективно улучшено алгоритмом Ньютона, как это происходило в оригинальной реализации для обратного квадратного корня. Я не изучал эту тему более глубоко и это, вероятно, будет темой отдельной публикации (наверное, не моей).
Выражение выше содержит новую “магическую константу” 0x3f7a3bea. В некотором смысле (из-за своей универсальности) она даже “более магическая”, чем константа в оригинальном коде. Давайте назовём её С и посмотрим на неё чуть внимательнее.
Давайте проверим работу нашей формулы для случая, когда p = 0. Как вы помните из курса математики, любое число в степени 0 равно единице. Что же будет с нашей формулой? Всё очень просто – умножение на 0 уничтожит второе слагаемое и у нас останется:
i = 0x3f7a3bea;
Что, действительно является константой и, будучи переведённой во float-формат, даст нам 0.977477 – т.е. “почти 1”. Поскольку мы имеем дело с аппроксимациями – это неплохое приближение. Кроме того, это говорит нам кое-что ещё. Наша константа С имеет совершенно не случайное значение. Это единица в формате чисел с плавающей запятой (ну или “почти единица”)
Это интересно. Давайте взглянем поближе:
Целочисленное представление C это:
![image](https://habrastorage.org/getpro/habr/post_images/8f8/04f/625/8f804f625d5a14378881ce21f12d1129.png)
Это почти, но всё же не совсем, форма числа с плавающей запятой. Единственная проблема – мы вычитаем вторую часть выражения, а должны были бы её прибавлять. Но это можно исправить:
![image](https://habrastorage.org/getpro/habr/post_images/54d/dc1/4e4/54ddc14e43de0edefb16a7e8fcef704f.png)
Вот теперь это выглядит в точности как целочисленное представление числа с плавающей запятой. Чтобы определить, какого именно числа, мы вычислим экспоненту и мантиссу, а потом уже и само число С. Вот экспонента:
![image](https://habrastorage.org/getpro/habr/post_images/181/499/46d/18149946d0b62fd204ecf880a65c6744.png)
А вот мантисса:
![image](https://habrastorage.org/getpro/habr/post_images/75f/f7d/c02/75ff7dc02e9cea442e4fa32953a9852b.png)
А, значит, значение самого числа будет равно:
![image](https://habrastorage.org/getpro/habr/post_images/1cb/c9a/34a/1cbc9a34acef9e94c7e349ffb35ece89.png)
И правда, если поделить наше σ (а оно было равно 0.0450465) на 2 и отнять результат от единицы, то мы получим уже известное нам 0.97747675, то самое, которое “почти 1”. Это позволяет нам посмотреть на С с другой стороны и вычислять его на рантайме:
float sigma = 0.0450465;
float c_sigma = 1 - (0.5f * sigma);
int C_sigma = *(*int)&c_sigma;
Учтите, что для фиксированного σ все эти числа будут константами и компилятор сможет рассчитать их на этапе компиляции. Результатом будет 0x3f7a3beb, что не совсем 0x3f7a3bea из расчетов выше, но отличается от него всего на 1 бит (наименее значимый). Получить из этого числа оригинальную константу из кода (и заголовка данной статьи) можно, умножив результат вычислений ещё на 1.5.
Со всеми этими выкладками мы приблизились к пониманию того, что в коде в начале статьи нет никакой “магии”, “фокусов”, “интуиции”, “подбора” и прочих грязных хаков, а есть лишь чистая математика, во всей своей первозданной красоте. Для меня главным выводом из этой истории стала новость о том, что преобразование float в int и наоборот путем переинтерпретации одного и того же набора бит – это не ошибка программиста и не “хак”, а вполне себе иногда разумная операция. Слегка, конечно, экзотическая, но зато очень быстрая и дающая практически полезные результаты. И, мне кажется, для этой операции могут быть найдены и другие применения – будем ждать.