БЫСТРОЕ ВЫЧИСЛЕНИЕ ОБРАТНОГО КВАДРАТНОГО КОРНЯ В 3D ПРОГРАММИРОВАНИИ

Одна из основных часто используемых операций в 3D программирование - это расчет обратного квадратного корня. При работе с 3D графикой используются нормали поверхностей (вектор длиной в единицу), которые нужны для описания освещения, отражения и т.д. Таких векторов нужно очень много, за секунды должно вычисляться миллионы нормированных векторов. Для нормализации вектора необходимо каждую координату вектора разделить на текущую длину вектора. Или если перефразировать каждую координату умножить на величину (Формула 1).

Формула 1. Обратная длинна вектора

 – расчёт этой части довольно прост, и все процессоры быстро считают это. А расчёт квадратного корня и операция деления, уже занимает больше времени на выполнение. Алгоритм представлен на С языке (Рисунок 1) и сейчас разъясним особенности этого алгоритма.

Рисунок 1. Код программы

Алгоритм состоит из четырех основных шагов. Первым делом вычисляем половину исходного числа, для метода Ньютона и сохраняем (float xhalf = 0.5f  * x;). Затем входное число которое пришло нам в формате float и интерпретируем его биты как значение типа integer (int i = *(int*)&x;). Далее над полученным целом числом производится довольно странная с первого взгляда операция, которая сдвигает число на один бит вправо и вычитает константу 5F3759DF  которую чуть позже поймем. Эта операция работает очень быстро и дает некую аппроксимацию требуемого результата (i = 0x5f3759df – (i >> 1);). Но данное число еще не является, собственно, результатом. Это лишь целое число, биты которого представляют некоторое число с плавающей запятой. Нам необходимо произвести обратное преобразование из integer в float (x = *(float*)&i;). И в итоги выполняется одна итерация метода Ньютона для улучшения аппроксимации (x = x * (1.5f – (xhalf*x*x));). В итоги мы имеем отличный алгоритм аппроксимации обратного квадратного корня.

Данный метод не подходит для точных вычислений, но для 3D программирования является отличным способом. Значения получаются довольно точные, но с погрешностью (+0% -0,17%), что для целей 3D графики вполне подходит.

Для выяснения странной операции нам необходимо вспомнить как хранятся числа с плавающей запятой. Рассматривая число с плавающей запятой как набор битов, экспонента и мантисса могут восприниматься как два положительных целых числа, обозначим их соответственно E и M. Интерпретируя биты числа с плавающей запятой, мы будем рассматривать мантиссу как число от 0 до 1. Экспоненту обозначим как знаковое число, вычтем смещение и получим число в диапазоне от -127 до 128.  И обозначим интерпретации этих значений как  e и m (Формула 2).

Формула 2. Интерпретация экспоненты и мантиссы

Где  и . Имея некоторые значения e и m, можно получить само число (Формула 3).

Формула 3. Представление числа

Нам приходит некоторое число X и требуется рассчитать его обратный квадратный корень (Формула 4).

Формула 4. Обратный квадратный корень

Возьмем логарифм с обоих сторон и представим X и Y как числа с плавающей запятой (Формула 5).

Формула 5. Представление формулы с плавающей запятой

Так как мантисса у нас в диапазоне от 0 до 1, то функция  достаточна близка к прямой линии и можно преобразовать в (Формула 6), где  некая константа.

Формула 6. Приближенная функция логарифма

Теперь перейдем от интерпретаций мантиссы и экспоненты к целочисленным представлениям, и выполнив несколько тривиальных преобразований, мы получим нечто довольно похожее на нашу строчку из кода (Формула 7).

Формула 7. Преобразованная формула без логарифмов

Ну в итоге если присмотреться то можно увидеть представление чисел с плавающей запятой из формулы что вывели раньше (Формула 8).

Формула 8. Конечное преобразование.

Теперь эта формула похожа на нашу строчку из кода (i = 0x5f3759df – (i >> 1);), но нам теперь надо найти константу. Мы уже знаем значения B и L, но нужно подобрать значение . Представим , так как это довольно не плохое приближение, и используя его мы получим (Формула 9).

Формула 9. Расчёт константы

В итоги это значение представим как HEX и получим 1597463007 = 5F3759DF. В итоге наг алгоритм рассчитывает очень быстро обратный квадратный корень.

Список литературы:

1. Амосов, А. А. Вычислительные методы: учебное пособие / А. А. Амосов, Ю. А. Дубинский, Н. В. Копченова. - 4-е изд., стер. - СПб. : Лань, 2014.
2. Mark Giambruno. 3D Grpahics & Animation Second Edition: New Riders Press, 2002