stringtranslate.com

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

В расчетах освещения и отражения, как в видеоигре OpenArena , используется быстрый код обратного квадратного корня для вычисления углов падения и отражения.

Быстрый обратный квадратный корень , иногда называемый Fast InvSqrt() или шестнадцатеричной константой 0x5F3759DF , представляет собой алгоритм, который оценивает обратную (или мультипликативную обратную) величину квадратного корня 32-битного числа с плавающей запятой в стандарте IEEE 754 с плавающей запятой. -точечный формат . Алгоритм наиболее известен благодаря своей реализации в 1999 году в Quake III Arena , видеоигре- шутере от первого лица, в значительной степени основанной на 3D-графике . С последующим развитием аппаратного обеспечения, особенно с инструкцией x86 SSE , этот алгоритм, как правило, не является лучшим выбором для современных компьютеров [1] , хотя он остается интересным историческим примером. [2] rsqrtss

Алгоритм принимает на вход 32-битное число с плавающей запятой и сохраняет уменьшенное пополам значение для дальнейшего использования. Затем, рассматривая биты, представляющие число с плавающей запятой, как 32-битное целое число, выполняется логический сдвиг вправо на один бит и результат вычитается из числа 0x 5F3759DF , которое является представлением с плавающей запятой аппроксимации . [3] Это приводит к первому приближению обратного квадратного корня входных данных. Снова рассматривая биты как числа с плавающей запятой, он запускает одну итерацию метода Ньютона , давая более точное приближение.

История

Уильям Кахан и К.С. Нг из Беркли в мае 1986 года написали неопубликованную статью, описывающую, как вычислить квадратный корень с использованием методов битовой игры с последующими итерациями Ньютона. [4] В конце 1980-х годов Клив Молер из Ardent Computer узнал об этой технике [5] и передал ее своему коллеге Грегу Уолшу. Грег Уолш разработал ныне известный константный и быстрый алгоритм обратного квадратного корня. Гэри Таролли консультировал Kubota, компанию, которая в то время финансировала Ardent, и, вероятно, представил алгоритм в 3dfx Interactive примерно в 1994 году. [6] [7]

Джим Блинн продемонстрировал простую аппроксимацию обратного квадратного корня в колонке журнала IEEE Computer Graphics and Applications в 1997 году . [8] Реверс-инжиниринг других современных 3D-видеоигр выявил вариацию алгоритма в игре Activision Interstate '76 1997 года . [9]

Quake III Arena , видеоигра- шутер от первого лица , была выпущена в 1999 году компанией id Software и использовала этот алгоритм. Брайан Хук, возможно, перенес этот алгоритм из 3dfx в id Software. [6] Обсуждение кода появилось на китайском форуме разработчиков CSDN в 2000 году, [10] а Usenet и форум gamedev.net широко распространили код в 2002 и 2003 годах. [11] Возникли предположения относительно того, кто написал алгоритм и как была получена константа; некоторые предполагали, что это Джон Кармак . [7] Полный исходный код Quake III был представлен на QuakeCon 2005 , но ответов на него не было. Ответ на вопрос об авторстве был получен в 2006 году, когда Грег Уолш связался с Beyond3D, поскольку их предположение приобрело популярность на Slashdot. [6]

В 2007 году алгоритм был реализован в некоторых специализированных аппаратных вершинных шейдерах с использованием программируемых пользователем вентильных матриц (FPGA). [12] [13]

Мотивация

Нормали поверхности широко используются в расчетах освещения и затенения, что требует расчета норм для векторов. Здесь показано поле векторов, нормальных к поверхности.
Двумерный пример использования нормали для нахождения угла отражения по углу падения; в данном случае — на свете, отражающемся от изогнутого зеркала. Быстрый обратный квадратный корень используется для обобщения этого расчета на трехмерное пространство.

Обратный квадратный корень из числа с плавающей запятой используется в цифровой обработке сигналов для нормализации вектора, масштабирования его до длины 1 для получения единичного вектора . [14] Например, программы компьютерной графики используют обратные квадратные корни для вычисления углов падения и отражения для освещения и затенения . Программы трехмерной графики должны выполнять миллионы таких вычислений каждую секунду, чтобы имитировать освещение. Когда код был разработан в начале 1990-х годов, большая часть вычислительной мощности с плавающей запятой отставала от скорости обработки целых чисел. [7] Это было проблематично для программ 3D-графики до появления специализированного оборудования для обработки преобразований и освещения . Вычисление квадратных корней обычно зависит от множества операций деления, которые для чисел с плавающей запятой требуют больших вычислительных затрат . Быстрый обратный квадрат дает хорошее приближение всего за один шаг деления.

Длина вектора определяется путем вычисления его евклидовой нормы : квадратного корня из суммы квадратов компонент вектора . Когда каждый компонент вектора делится на эту длину, новый вектор будет единичным вектором , указывающим в том же направлении. В программе трехмерной графики все векторы находятся в трехмерном пространстве , поэтому вектор .

– евклидова норма вектора.

– нормированный (единичный) вектор. Подставив , уравнение можно также записать как:

который связывает единичный вектор с обратным квадратным корнем компонентов расстояния. Обратный квадратный корень можно использовать для вычислений, поскольку это уравнение эквивалентно

где член дроби представляет собой обратный квадратный корень из .

В то время деление чисел с плавающей запятой было, как правило, дорогостоящим по сравнению с умножением; быстрый алгоритм обратного квадратного корня обошел этап деления, что дало ему преимущество в производительности.

Обзор кода

Следующий код представляет собой быструю реализацию обратного квадратного корня из Quake III Arena , лишенную директив препроцессора C , но включающую точный исходный текст комментария: [15]

float q_rsqrt ( число с плавающей запятой ) { long i ; плавать x2 , y ; const float Threehalfs = 1.5F ;             х2 = число * 0,5F ; у = число ; я = * ( длинный * ) & y ; // злой взлом битового уровня с плавающей запятой i = 0x5f3759df - ( i >> 1 ); // что за херня? y = * ( float * ) & я ; y = y * ( три половинки - ( x2 * y * y ) ); // 1-я итерация // y = y * (threehalfs - (x2 * y * y)); // 2-я итерация, это можно удалить                                                    вернуть y ; } 

В то время общий метод вычисления обратного квадратного корня заключался в вычислении аппроксимации для , а затем пересмотре этой аппроксимации другим методом до тех пор, пока она не попадала в допустимый диапазон ошибок фактического результата. Распространенные программные методы в начале 1990-х годов основывались на справочной таблице . [16] Ключом к быстрому обратному квадратному корню было непосредственное вычисление аппроксимации с использованием структуры чисел с плавающей запятой, что оказалось быстрее, чем поиск в таблице. Алгоритм оказался примерно в четыре раза быстрее, чем вычисление квадратного корня другим методом и вычисление обратной величины посредством деления с плавающей запятой. [17] Алгоритм был разработан с учетом спецификации 32-битных чисел с плавающей запятой IEEE 754-1985 , но исследование Криса Ломонта показало, что он может быть реализован и в других спецификациях с плавающей запятой. [18]

Преимущества в скорости, обеспечиваемые быстрым трюком с обратным квадратным корнем, заключаются в том, что 32-битное слово с плавающей запятой [примечание 1] рассматривается как целое число , а затем вычитается из « магической » константы 0x 5F3759DF . [7] [19] [20] [21] Это вычитание целого числа и сдвиг битов приводит к образованию битового шаблона, который, если его переопределить как число с плавающей запятой, является грубым приближением обратного квадратного корня из числа. Для достижения некоторой точности выполняется одна итерация метода Ньютона, и код готов. Алгоритм генерирует достаточно точные результаты, используя уникальное первое приближение метода Ньютона ; однако это намного медленнее и менее точно, чем использование инструкции SSE на процессорах x86, также выпущенных в 1999 году. [1] [22]rsqrtss

Рабочий пример

Например, число можно использовать для расчета . Первые шаги алгоритма показаны ниже:

0011_1110_0010_0000_0000_0000_0000_0000 Битовая комбинация x и i0001_1111_0001_0000_0000_0000_0000_0000 Сдвиг вправо на одну позицию: (i >> 1)0101_1111_0011_0111_0101_1001_1101_1111 Магическое число 0x5F3759DF0100_0000_0010_0111_0101_1001_1101_1111 Результат 0x5F3759DF - (i >> 1)

Интерпретация как 32-битное представление IEEE:

0_01111100_0100000000000000000000 1,25 × 2 −3
0_00111110_00100000000000000000000 1,125 × 2 −65
0_10111110_01101110101100111011 111 1,432430... × 2 63
0_10000000_01001110101100111011111 1,307430... × 2 1

Переинтерпретация этого последнего битового шаблона как числа с плавающей запятой дает аппроксимацию с погрешностью около 3,4%. После одной-единственной итерации метода Ньютона окончательный результат составляет погрешность всего 0,17%.

Как избежать неопределенного поведения

Согласно стандарту C, переинтерпретация значения с плавающей запятой как целого числа путем приведения и разыменования указателя на него недопустима ( неопределенное поведение ). Другой способ — поместить значение с плавающей запятой в анонимное объединение , содержащее дополнительный 32-битный элемент целого числа без знака, и доступ к этому целому числу обеспечивает представление на уровне битов содержимого значения с плавающей запятой. Однако игра слов типов через объединение также является неопределенным поведением в C++.

# включаем <stdint.h> // uint32_t  float q_rsqrt ( число с плавающей запятой ) { union { float f ; uint32_t я ; } Конв = { . е = число }; конв . я = 0x5f3759df - ( услов . я >> 1 ); конв . f *= 1,5F - ( число * 0,5F * услов . f * услов . f ); возврат конв . ж ; }                                    

В современном C++ обычным методом реализации приведения этой функции является использование C++20 std::bit_cast. Это также позволяет функции работать в constexprконтексте:

# включаем <бит> # включаем <лимиты> # включаем <cstdint>      constexpr float Q_rsqrt ( число с плавающей запятой ) noException { static_assert ( std :: numeric_limits <float> :: is_iec559 ); // (включить только в IEEE 754)       float const y = std :: bit_cast <float> ( 0x5f3759df - ( std :: bit_cast < std :: uint32_t > ( число ) >> 1 ) ) ; вернуть y * ( 1,5f - ( число * 0,5f * y * y )); }                     

Алгоритм

Алгоритм выполняет вычисления , выполняя следующие шаги:

  1. Псевдоним аргумента целого числа как способ вычисления приближения двоичного логарифма.
  2. Используйте это приближение, чтобы вычислить приближение
  3. Псевдоним обратного значения с плавающей запятой, как способ вычисления аппроксимации экспоненты по основанию 2.
  4. Уточните приближение, используя одну итерацию метода Ньютона.

Представление с плавающей запятой

Поскольку этот алгоритм в значительной степени опирается на представление чисел с плавающей запятой одинарной точности на уровне битов, здесь представлен краткий обзор этого представления. Чтобы закодировать ненулевое действительное число как число с плавающей запятой одинарной точности, первым шагом является запись в виде нормализованного двоичного числа: [23]

где показатель степени является целым числом и является двоичным представлением мантиссы . Поскольку единственный бит перед точкой в ​​мантиссе всегда равен 1, его не нужно сохранять. Уравнение можно переписать как:

где значит , поэтому . Из этой формы вычисляются три целых числа без знака: [24]

Затем эти поля упаковываются слева направо в 32-битный контейнер. [25]

В качестве примера снова рассмотрим число . Нормализация урожайности:

и, таким образом, три целочисленных поля без знака:

эти поля упакованы, как показано на рисунке ниже:

Псевдоним целого числа как приближенный логарифм

Если бы расчеты нужно было производить без компьютера или калькулятора, была бы полезна таблица логарифмов вместе с тождеством , действительным для любого основания . Быстрый обратный квадратный корень основан на этом тождестве, а также на том факте, что псевдоним float32 для целого числа дает грубое приближение его логарифма. Вот как:

Если — положительное нормальное число :

затем

и поскольку , логарифм в правой части можно аппроксимировать формулой [26]

где – свободный параметр, используемый для настройки аппроксимации. Например, дает точные результаты на обоих концах интервала, в то же время дает оптимальное приближение (лучшее в смысле равномерной нормы ошибки). Однако это значение не используется алгоритмом, поскольку не учитывает последующие шаги.

Целое число, связанное с числом с плавающей запятой (синим цветом), по сравнению с масштабированным и сдвинутым логарифмом (серым цветом).

Таким образом, существует приближение

Интерпретация битового шаблона с плавающей запятой как целого числа дает [примечание 4]

Тогда оказывается, что это масштабированная и сдвинутая кусочно-линейная аппроксимация , как показано на рисунке справа. Другими словами, аппроксимируется

Первое приближение результата

Расчет основан на тождестве

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

Таким образом, приближение равно :

который записан в коде как

я = 0x5f3759df - ( я >> 1 );        

Первое слагаемое выше — это магическое число

из чего можно сделать вывод . Второе слагаемое вычисляется путем сдвига битов на одну позицию вправо. [27]

метод Ньютона

Относительная ошибка между прямым вычислением и быстрым обратным квадратным корнем при выполнении 0, 1, 2, 3 и 4 итераций метода поиска корня Ньютона. Обратите внимание, что используется двойная точность, и наименьшая представимая разница между двумя числами двойной точности достигается после выполнения 4 итераций.

С обратным квадратным корнем . Аппроксимация, полученная на предыдущих этапах, может быть уточнена с помощью метода поиска корня — метода, который находит ноль функции . Алгоритм использует метод Ньютона : если существует приближение для , то лучшее приближение можно вычислить, взяв , где – производная от at . [28] В связи с вышеизложенным ,

где и .

Обработка как числа с плавающей запятой эквивалентнаy = y*(threehalfs - x/2*y*y);

Повторяя этот шаг, используя выходные данные функции ( ) в качестве входных данных следующей итерации, алгоритм приводит к сходимости к обратному квадратному корню. [29] Для целей движка Quake III использовалась только одна итерация. Вторая итерация осталась в коде, но была закомментирована . [21]

Точность

Как отмечалось выше, приближение очень точное. На единственном графике справа показана ошибка функции (то есть ошибка аппроксимации после того, как она была улучшена за счет выполнения одной итерации метода Ньютона) для входных данных, начинающихся с 0,01, где стандартная библиотека дает в результате 10,0. , а InvSqrt() дает 9,982522, что составляет относительную разницу 0,0017478, или 0,175% от истинного значения, равного 10. С этого момента абсолютная ошибка только падает, а относительная ошибка остается в тех же пределах во всех порядках величины.

Последующие улучшения

Магическое число

Точно неизвестно, как было определено точное значение магического числа. Крис Ломонт разработал функцию, позволяющую минимизировать ошибку аппроксимации , выбирая магическое число в диапазоне. Сначала он вычислил оптимальную константу для шага линейной аппроксимации как 0x5F37642F , близкую к 0x5F3759DF , но эта новая константа дала немного меньшую точность после одной итерации метода Ньютона. [30] Затем Ломонт искал постоянный оптимум даже после одной и двух итераций Ньютона и нашел 0x5F375A86 , что более точно, чем оригинал на каждом этапе итерации. [30] В заключение он задал вопрос, было ли точное значение исходной константы выбрано путем вывода или методом проб и ошибок . [31] Ломонт сказал, что магическое число для 64-битного типа размера IEEE754 double — это 0x5FE6EC85E7DE30DA , но позже Мэтью Робертсон показал, что оно равно именно 0x5FE6EB50C7B537A9 . [32]

Ян Кадлец уменьшил относительную ошибку еще в 2,7 раза, корректируя константы и в единственной итерации метода Ньютона, [33] придя после исчерпывающего поиска к

конв . я = 0x5F1FFFF9 - ( услов . я >> 1 ); конв . f *= 0,703952253f * ( 2,38924456f - x * усл . f * усл . f ); возврат конв . ж ;                     

Полный математический анализ для определения магического числа теперь доступен для чисел с плавающей запятой одинарной точности. [34] [35]

Нулевой результат

Промежуточным по сравнению с использованием одной и двух итераций метода Ньютона с точки зрения скорости и точности является одна итерация метода Галлея . В этом случае метод Галлея эквивалентен применению метода Ньютона с исходной формулой . Затем этап обновления

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

Устаревание

Последующие дополнения производителей оборудования сделали этот алгоритм по большей части излишним. Например, для x86 компания Intel представила инструкцию SSErsqrtss в 1999 году. В тесте Intel Core 2 2009 года эта инструкция занимала 0,85 нс на одно число с плавающей запятой по сравнению с 3,54 нс для быстрого алгоритма обратного квадратного корня и имела меньшую ошибку. [1]

Некоторые недорогие встроенные системы не имеют специализированных инструкций по извлечению квадратного корня. Однако производители этих систем обычно предоставляют тригонометрические и другие математические библиотеки, основанные на таких алгоритмах, как CORDIC . Большинство приложений в этой области не чувствительны к производительности, и специализированный быстрый алгоритм обратного квадратного корня не требуется. [ нужна цитата ]

Смотрите также

Примечания

  1. ^ Использование этого типа longснижает переносимость этого кода в современных системах. Для правильного выполнения кода размер sizeof(long)должен составлять 4 байта, в противном случае могут возникнуть отрицательные выходные данные. Во многих современных 64-битных системах sizeof(long)он составляет 8 байт. Более портативная замена — int32_t.
  2. ^ должен находиться в диапазоне , чтобы его можно было представить как обычное число .
  3. ^ Единственные действительные числа, которые могут быть представлены точно как с плавающей запятой, - это те, для которых является целым числом. Остальные числа можно представить только приблизительно, округлив их до ближайшего точно представимого числа.
  4. ^ Поскольку положительно, .

Рекомендации

  1. ^ abc Раскин, Элан (16 октября 2009 г.). «Квадратный корень времени». Требуется некоторая сборка . Архивировано из оригинала 08 февраля 2021 г. Проверено 7 мая 2015 г.
  2. ^ фейлипу. «z88dk — это набор инструментов разработки программного обеспечения, предназначенный для компьютеров 8080 и z80». Гитхаб .
  3. ^ Мунафо, Роберт. «Заметные свойства конкретных чисел». mrob.com . Архивировано из оригинала 16 ноября 2018 года.
  4. ^ «Реализация sqrt в fdlibm — см. обсуждение У. Кахана и К.К. Нг в комментариях в нижней половине этого кода».
  5. Молер, Клив (19 июня 2012 г.). «Симплектическая космическая война». MATLAB Central — Уголок Клива . МАТЛАБ . Проверено 21 июля 2014 г.
  6. ^ abc Соммефельдт, Рис (19 декабря 2006 г.). «Происхождение Fast InvSqrt() в Quake3 — Часть вторая». За пределами 3D . Проверено 19 апреля 2008 г.
  7. ^ abcd Соммефельдт, Рис (29 ноября 2006 г.). «Происхождение Fast InvSqrt() в Quake3». За пределами 3D . Проверено 12 февраля 2009 г.
  8. ^ Блинн 1997, стр. 80–84.
  9. Пилар, Шейн (1 июня 2021 г.). «Быстрый обратный квадратный корень... в 1997 году?!».
  10. ^ «Обсуждение на CSDN» . Архивировано из оригинала 2 июля 2015 г.
  11. ^ Ломонт 2003, с. 1-2.
  12. ^ Зафар, Саад; Адапа, Равитея (январь 2014 г.). «Проектирование аппаратной архитектуры и отображение алгоритма быстрого обратного квадратного корня». Международная конференция по достижениям в электротехнике (ICAEE) , 2014 г. стр. 1–4. doi : 10.1109/ICAEE.2014.6838433. ISBN 978-1-4799-3543-7. S2CID  2005623.
  13. ^ Миддендорф 2007, стр. 155–164.
  14. ^ Блинн 2003, с. 130.
  15. ^ "quake3-1.32b/code/game/q_math.c". Квейк III Арена . программное обеспечение id . Архивировано из оригинала 2 августа 2020 г. Проверено 21 января 2017 г.
  16. ^ Эберли 2001, с. 504.
  17. ^ Ломонт 2003, с. 1.
  18. ^ Ломонт 2003.
  19. ^ Ломонт 2003, с. 3.
  20. ^ МакЭнири 2007, с. 2, 16.
  21. ^ аб Эберли 2001, с. 2.
  22. ^ Туман, Агнер. «Списки задержек инструкций, пропускной способности и сбоев в микрооперациях для процессоров Intel, AMD и VIA» (PDF) . Проверено 8 сентября 2017 г.
  23. ^ Голдберг 1991, с. 7.
  24. ^ Голдберг 1991, стр. 15–20.
  25. ^ Голдберг 1991, с. 16.
  26. ^ МакЭнири 2007, с. 3.
  27. ^ Хеннесси и Паттерсон 1998, с. 305.
  28. ^ Харди 1908, с. 323.
  29. ^ МакЭнири 2007, с. 6.
  30. ^ аб Ломонт 2003, с. 10.
  31. ^ Ломонт 2003, стр. 10–11.
  32. ^ Мэтью Робертсон (24 апреля 2012 г.). «Краткая история InvSqrt» (PDF) . УНБСЖ.
  33. ^ Кадлец, Январь (2010). "Řrřlog::Улучшение быстрого извлечения квадратного корня" (личный блог) . Проверено 14 декабря 2020 г.
  34. ^ Мороз и др. 2018.
  35. Мюллер, Жан-Мишель (декабрь 2020 г.). «Элементарные функции и приближенные вычисления». Труды IEEE . 108 (12): 2146. doi :10.1109/JPROC.2020.2991885. ISSN  0018-9219. S2CID  219047769.

Библиография

дальнейшее чтение

Внешние ссылки