Быстрый обратный квадратный корень: различия между версиями
[непроверенная версия] | [отпатрулированная версия] |
Содержимое удалено Содержимое добавлено
Mercury (обсуждение | вклад) |
РобоСтася (обсуждение | вклад) м чистка управляющих символов Юникода |
||
(не показаны 44 промежуточные версии 2 участников) | |||
Строка 1:
[[Файл:OpenArena-Rocket.jpg|thumb|300px|right|При расчёте освещения ''[[OpenArena]]'' (свободный порт ''[[Quake III: Arena]]'') вычисляет [[угол падения|углы падения]] и отражения через быстрый обратный квадратный корень. Обратите внимание на кожух оружия — при очень низкой детализации (8 четырёхугольников) игра делает вид, что он криволинейный.]]
'''Бы́стрый обра́тный квадра́тный ко́рень''' (также ''быстрый InvSqrt()'' или ''0x5F3759DF'' по используемой [[Магическое число (программирование)|«магической» константе]]
== Алгоритм ==
Алгоритм принимает 32-битное число с плавающей запятой ([[число одинарной точности|одинарной точности]] в формате [[IEEE 754]]) в качестве исходных данных и производит над ним следующие операции:
# Трактуя 32-битное дробное число как целое, провести операцию {{nobr|1=y<sub>0</sub> = 5F3759DF[[Шестнадцатеричная система счисления|<sub>16</sub>]] − (x >> 1)}}, где >> — [[битовый сдвиг]] вправо. Результат снова трактуется как 32-битное дробное число.
# Для уточнения можно провести одну итерацию [[Метод Ньютона|метода Ньютона]]: {{nobr|1 =
Реализация из Quake III<ref name="hummus" />:
<syntaxhighlight lang="c" line="1">
Строка 17:
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; //
i = 0x5f3759df - ( i >> 1 ); //
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); //
// y = y * ( threehalfs - ( x2 * y * y ) ); //
return y;
}
</syntaxhighlight>
Эта реализация считает, что {{cpp|float}} по длине равен {{cpp|long}}, и использует для преобразования указатели (может ошибочно сработать оптимизация «если изменился {{cpp|float}}, ни один {{cpp|long}} не менялся»; на GCC при компиляции в «выпуск» срабатывает предупреждение). По комментариям видно, что [[Джон Кармак]], выкладывая игру в открытый доступ, не понял, что там делается
Корректная по меркам современного Си реализация, с учётом возможных оптимизаций и [[кроссплатформенность|кроссплатформенности]]:
Строка 40:
float f;
uint32_t i;
} conv = {number}; //
conv.i = 0x5f3759df - ( conv.i >> 1 );
conv.f *= threehalfs - x2 * conv.f * conv.f;
Строка 68:
Саму идею приближения дробного числа целым для вычисления <math>\sqrt x</math> придумали [[Кэхэн, Уильям Мортон|Уильям Кэхэн]] и К. Ын в 1986<ref>{{Cite web |url=http://www.netlib.org/fdlibm/e_sqrt.c |title=Источник |access-date=2023-04-09 |archive-date=2023-04-13 |archive-url=https://web.archive.org/web/20230413082942/https://netlib.org/fdlibm/e_sqrt.c |deadlink=no }}</ref>. До этой идеи добрались Грег Уолш, [[Клив Моулер]] и Гэри Таролли, работавшие тогда в [[Stardent Inc.|Ardent Computer]]<ref>{{Cite web |url=https://www.beyond3d.com/content/articles/15/ |title=Beyond3D — Origin of Quake3’s Fast InvSqrt() — Part Two<!-- Заголовок добавлен ботом --> |access-date=2019-08-25 |archive-date=2019-08-25 |archive-url=https://web.archive.org/web/20190825105752/https://www.beyond3d.com/content/articles/15/ |deadlink=no }}</ref><ref>{{Cite web |url=https://blogs.mathworks.com/cleve/2012/06/19/symplectic-spacewar/#comment-13 |title=Symplectic Spacewar » Cleve’s Corner: Cleve Moler on Mathematics and Computing - MATLAB & Simulink<!-- Заголовок добавлен ботом --> |access-date=2023-04-09 |archive-date=2023-04-09 |archive-url=https://web.archive.org/web/20230409065120/https://blogs.mathworks.com/cleve/2012/06/19/symplectic-spacewar/#comment-13 |deadlink=no }}</ref>. Грегу Уолшу и приписывается знаменитая константа 0x5F3759DF.
Таролли, перешедший в [[3dfx]], принёс алгоритм туда, где он и применялся для расчёта углов падения и отражения света в трёхмерной графике. Джим Блинн, специалист по 3D-графике, переизобрёл метод в 1997 году с более простой константой 1,5<ref name="blinn">{{Cite web |url=https://ieeexplore.ieee.org/document/595279/ |title=Floating-point tricks {{!}} IEEE Journals & Magazine {{!}} IEEE Xplore<!-- Заголовок добавлен ботом --> |access-date=2022-08-17 |archive-date=2022-08-17 |archive-url=https://web.archive.org/web/20220817144947/https://ieeexplore.ieee.org/document/595279/ |deadlink=no }}. [https://www.yumpu.com/en/document/view/6104114/floating-point-tricks-ieee-computer-graphics-and-applications Копия на Yumpu] {{Wayback|url=https://www.yumpu.com/en/document/view/6104114/floating-point-tricks-ieee-computer-graphics-and-applications |date=20230409061404 }}</ref>. Более сложный табличный метод, который считает до 4 знаков (0,01 %), найден при дизассемблировании игры [[Interstate ’76]] (1997)<ref name="i76">{{Cite web |url=https://inbetweennames.net/blog/2021-05-06-i76rsqrt/ |title=Fast reciprocal square root… in 1997?! — Shane Peelar’s Blog<!-- Заголовок добавлен ботом --> |access-date=2022-08-17 |archive-date=2022-10-11 |archive-url=https://web.archive.org/web/20221011200325/https://inbetweennames.net/blog/2021-05-06-i76rsqrt/ |deadlink=no }}</ref>.
Однако данный метод не появлялся на общедоступных форумах, таких как [[Usenet]], до 2002—2003-х годов. Метод обнаружили в [[Quake III: Arena]], опубликованном в 2005, и приписали авторство [[Кармак, Джон|Джону Кармаку]]. Тот предположил, что его в [[id Software]] принёс [[Майкл Абраш]], специалист по графике, или Терье Матисен, специалист по ассемблеру; другие ссылались на Брайана Хука, выходца из 3dfx<ref>{{Cite web |url=https://www.beyond3d.com/content/articles/8/ |title=Beyond3D — Origin of Quake3’s Fast InvSqrt()<!-- Заголовок добавлен ботом --> |access-date=2019-10-04 |archive-date=2017-04-10 |archive-url=https://web.archive.org/web/20170410090337/https://www.beyond3d.com/content/articles/8/ |deadlink=no }}</ref>. Изучение вопроса показало, что код имел более глубокие корни как в аппаратной, так и в программной сферах компьютерной графики. Исправления и изменения производились как Silicon Graphics, так и [[3dfx Interactive]], при этом самая ранняя известная версия написана Гэри Таролли для [[SGI Indigo]].
Строка 132:
|colspan=4 style="border: none; border-right: 1px solid gray; text-align: right"|0
|}
Имеем дело только с положительными числами (знаковый бит равен нулю), не [[Денормализованные числа|денормализованными]], не [[Машинная бесконечность|<big>∞</big>]] и не [[NaN]]. Такие числа в [[экспоненциальная запись|стандартном виде]] записываются как 1,mmmm<sub>2</sub>·2<sup>e</sup>. Часть 1,mmmm называется ''мантиссой'', ''e'' — ''порядком''. Головную единицу не хранят (''неявная единица''), так что величину 0,mmmm назовём ''явной частью'' мантиссы. Кроме того, у машинных дробных чисел ''смещённый порядок'': 2<sup>0</sup> записывается как 011.1111.1<sub>2</sub>{{efn|name="dots"|Здесь и далее точки — границы полубайтов, вертикальные линии — границы полей компьютерного дробного.}}.
На положительных числах [[биекция]] «дробное ↔ целое» (ниже обозначенная как <math>I_x</math>) непрерывна как [[кусочно-линейная функция]] и [[монотонная функция|монотонна]]. Отсюда сразу же можно заявить, что быстрый обратный корень, как комбинация непрерывных функций, непрерывен. А первая его часть — сдвиг-вычитание — к тому же монотонна и кусочно-линейна. Биекция сложна, но почти «бесплатна»: в зависимости от архитектуры процессора и [[соглашение о вызове|соглашений вызова]], нужно или ничего не делать, или переместить число из дробного регистра в целочисленный.
В примере выше целочисленное представление равняется 0x3E20.0000, и оно раскладывается так: знаковое поле 0, смещённый порядок{{efn|name="dots"}} 011.1110.0<sub>2</sub>=124, несмещённый 124−127=−3, мантисса (вместе с неявной единицей) 1,01<sub>2</sub>=1,25, и дробное значение {{nobr|1=1,25·2<sup>−3</sup>=0,15625}}.
Обозначим <math>m_x \in [0,1)</math> явную часть мантиссы числа <math>x</math>, <math>e_x \in \mathbb{Z}</math> — несмещённый порядок, <math>L=2^{23}</math> — разрядность мантиссы, <math>B=127</math> — смещение порядка. Число{{efn|Здесь и далее ≡ означает «равно по определению».}} <math>x \equiv 2^{e_x} (1 + m_x)</math> будет иметь целочисленное представление <math>I_x \equiv L(e_x + B + m_x)</math>. Можно выписать и обратное преобразование: <math>e_x = \left\lfloor \tfrac{I_x}{L} - B \right\rfloor</math>, <math>m_x = \left\{ \tfrac{I_x}{L} \right\}</math>, ведь первое целое, а второе от 0 до 1.
Строка 159:
Это и есть формула первого приближения.
=== Магическая константа {{math|Q}} ===
Магическая константа <math>Q \equiv \tfrac {3}{2} L (B - \sigma)</math> находится в пространстве компьютерных целых, но её дробное представление <math>I^{-1}(Q)</math> также важно для исследователей. Его ''несмещённый'' порядок
: <math>e_Q = \left\lfloor \tfrac{Q}{L} - B\right\rfloor = \left\lfloor \tfrac{3L}{2L} (B - \sigma) - B \right\rfloor = \left\lfloor \tfrac{B - 3\sigma}{2} \right\rfloor</math>.
Строка 167:
Смещение порядка {{math|''B''}} нечётное, и полная мантисса числа <math>I^{-1}(Q)</math> равняется
: <math>c \equiv 1+m_Q = 1 + \left\{ \tfrac{Q}{L} \right\} = 1 + 0{,}5 - \tfrac{3}{2} \sigma \in (1{,}37; 1{,}5</math>),
а в двоичной записи{{efn|name="dots"}} — 0|101.1111.0|01<sub>1</sub>… (1 — неявная единица; 0,5 пришли из порядка; маленькая единица соответствует диапазону [1,375; 1,5) и потому крайне вероятна, но не гарантирована нашими прикидочными расчётами.)
[[Файл:Fast inverse square root linear approximation.svg|thumb|upright=1.5|Первое (кусочно-линейное) приближение быстрого обратного квадратного корня {{nobr|({{math|''c''}} {{=}} 1,43)}}]]
Через константу {{math|c}} можно вычислить, чему равняется первое кусочно-линейное приближение<ref name="moroz" /> (в источнике используется не сама мантисса, а её явная часть <math>m_Q \equiv t = c - 1</math>):
* Для <math>x \in [0{,}5; \; c - 0{,}5)</math>: <math>y_{01} = -x + t + \tfrac{3}{2} = -x + c + \tfrac{1}{2}</math>;
* Для <math>x \in [c - 0{,}5; \; 1)</math>: <math>y_{02} = -\tfrac{1}{2}x + \tfrac{1}{2}t + \tfrac{5}{4} = -\tfrac{1}{2}x + \tfrac{1}{2}c + \tfrac{3}{4}</math>;
Строка 181:
[[Метод Ньютона]] даёт<ref name="moroz" /> <math>f(y)=\frac{1}{y^2}-x</math>, <math>f'(y)=-\frac{2}{y^3}</math>, и <math>y_{n+1} = y_{n} - \frac{f(y_n)}{f'(y_n)} = \frac{y_{n}(3-x y_n^2)}{2} = y_{n}(1{,}5-0{,}5 \, x y_n^2)</math>. Функция <math>f(y)</math> убывает и выпукла вниз, на таких функциях метод Ньютона подбирается к истинному значению слева — потому алгоритм всегда занижает ответ.
После одного шага метода Ньютона результат получается довольно точный ({{nobr|+
Метод Ньютона не гарантирует монотонности, но компьютерный перебор показывает, что монотонность всё-таки есть.
Строка 250:
== Мотивация ==
[[Файл:Lighting prism cylinder.svg|thumb|upright=1.5|Поле нормалей: {{nobr|а) для}} призмы (угловатый объект); {{nobr|б) для}} низкополигонального цилиндра (
«Прямое» наложение освещения на [[Полигональная сетка|трёхмерную модель]], даже высокополигональную, даже с учётом [[Закон Ламберта|закона Ламберта]] и других сложных формул отражения и рассеивания, сразу же выдаст полигональный вид — зритель увидит разницу в освещении по рёбрам многогранника<ref name="normal_map">{{Cite web |url=https://habr.com/ru/post/481480/ |title=Это норма: что такое карты нормалей и как они работают / Хабр<!-- Заголовок добавлен ботом --> |access-date=2022-07-04 |archive-date=2020-07-10 |archive-url=https://web.archive.org/web/20200710004038/https://habr.com/ru/post/481480/ |deadlink=no }}</ref>. Иногда так и нужно — если предмет действительно угловатый. А для криволинейных предметов поступают так: в трёхмерной программе указывают, острое ребро или сглаженное<ref name="normal_map" />. В зависимости от этого ещё при экспорте модели в игру по углам треугольников вычисляют нормаль [[единичный вектор|единичной длины]] к криволинейной поверхности. При анимации и поворотах игра преобразует эти нормали вместе с остальными трёхмерными данными; при наложении освещения — интерполирует по всему треугольнику и нормализует (доводит до единичной длины).
Чтобы нормализовать вектор, надо разделить все три его компонента на длину. Или, что лучше, умножить их на величину, обратную длине: <math>(x', y', z') = (x, y, z)\frac{1}{\sqrt{x^2 + y^2 + z^2}}</math>. За секунду должны вычисляться миллионы этих корней. До того как было создано специальное аппаратное обеспечение для обработки [[Transform, clipping, and lighting|трансформаций и освещения]], программное обеспечение вычислений могло быть медленным. В частности, в начале 1990-х, когда код был разработан, большинство вычислений с плавающей запятой отставало по производительности от операций с целыми числами.
Строка 262:
При желании можно перебалансировать погрешность, умножив коэффициенты 1,5 и 0,5 на 1,0009, чтобы метод давал симметрично ±0,09 % — так поступили<ref name="i76" /> в игре [[Interstate ’76]], которая также делает итерацию метода Ньютона.
Чех Ян Ка́длец [[двоичный поиск|двоичным поиском]], а затем перебором в окрестности найденного улучшил алгоритм<ref name="kadlec">{{Cite web |url=http://rrrola.wz.cz/inv_sqrt.html |title=Řrřlog :: Improving the fast inverse square root<!-- Заголовок добавлен ботом --> |access-date=2023-04-08 |archive-date=2023-04-08 |archive-url=https://web.archive.org/web/20230408222901/http://rrrola.wz.cz/inv_sqrt.html |deadlink=no }}</ref>. Его метод даёт в 1,3 раза меньшую симметричную погрешность — не ±0,09 %, а ±0,065.
|