$30 off During Our Annual Pro Sale. View Details »

0x5F3759DF

ykozw
October 15, 2024

 0x5F3759DF

Quake III Arenaのコードの中にあったことで注目された高速逆平方根アルゴリズム「Fast Inverse Square Root」を解説します。0x5F3759DFというマジックナンバーの正体と、この方法でなぜ逆平方根が計算できるのかを、ビット操作や数学的原理に触れながら詳しく説明します。
🔗: レイトレ合宿10(https://sites.google.com/view/rtcamp10)

ykozw

October 15, 2024
Tweet

More Decks by ykozw

Other Decks in Technology

Transcript

  1. 次の関数は何をする? float doSomething(float x) { const uint32_t t = 0x5f3759df

    - (bit_cast<uint32_t>(x) >> 1); const float y = bit_cast<float>(t); return y * (1.5f - (x * 0.5f * y * y)); } ←tをfloatと解釈する ← 何らかの変換? ↓??? 実は rsqrt(f x = 1 √𝑥 )の(変わった)実装!
  2. Fast inverse square root Quake III Arena(1999)にて実装されてい たことで注目を浴びる ベクトルの正規化で使われる つまりレイトレでも使える

    現代では`RSQRTSS`の方が高速で高精度 現代でも学ぶ価値がある(と思う) John Carmackが発案したわけではない "John Carmack E3 2006" by Rob Fahey is licensed under CC BY-SA 2.0.
  3. rsqrt float rsqrt(float x) { const uint32_t t = 0x5f3759df

    - (bit_cast<uint32_t>(x) >> 1); const float y = bit_cast<float>(t); return y * (1.5f - (x * 0.5f * y * y)); }
  4. y*(1.5f-(x*0.5f*y*y)) 次の求根問題の解が欲しい 1 √𝑥 − 𝑦 = 0 変形すると次のようになる 𝑓

    𝑦 = 1 𝑦2 − 𝑥 この𝑓 𝑦 = 0となるyを得たい。 微分すると次のようになる 𝑓′ 𝑦 = − 2 𝑦3 "Newton iteration" by Yves Daoust is licensed under Public Domain.
  5. y*(1.5f-(x*0.5f*y*y)) 新しい解 = 今の解 - f(今の解)/f‘(今の解) 𝑦 − 𝑓 𝑦

    𝑓′ 𝑦 = 𝑦 − 1 𝑦2 − 𝑥 − 2 𝑦3 = 𝒚( 𝟑 𝟐 − 𝒙𝒚𝟐 𝟐 ) "Newton iteration" by Yves Daoust is licensed under Public Domain.
  6. rsqrt float rsqrt(float x) { const uint32_t t = 0x5f3759df

    - (bit_cast<uint32_t>(x) >> 1); const float y = bit_cast<float>(t); return y * (1.5f - (x * 0.5f * y * y)); }
  7. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) 唐突にlog2 (𝑥)を取ることを考える log2 (𝑥) = log2 (2𝐸−127

    1 + 1 223 𝑀 ) = log2 2𝐸−127 + log2 1 + 1 223 𝑀 = 𝑬 − 𝟏𝟐𝟕 + 𝒍𝒐𝒈𝟐 𝟏 + 𝟏 𝟐𝟐𝟑 𝑴 "Float example" by Dave Barber is licensed under CC BY-SA 3.0.
  8. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) 一般に次が成り立つ(マクローリン展開) log2 1 + 𝑥 =𝑥 −

    𝑥2 2 + 𝑥3 3 − ⋯ 今回のx = 𝟏 𝟐𝟐𝟑 𝑴は[0,1] の範囲しかとらない 𝑙𝑜𝑔2 1 + 1 223 𝑀 大胆に次のように近似する log2 1 + 𝑥 ≈𝑥 + 𝜇 このとき、 𝜇を0.043にするとよく近似できる ※ 詳しくは後述
  9. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) log2 (𝑥) = 𝐸 − 127 +

    𝑙𝑜𝑔2 1 + 1 223 𝑀 ≈ 𝐸 − 127 + 1 223 𝑀 + 𝜇 = 1 223 223𝐸 + 𝑀 + 𝜇 − 127 "Float example" by Dave Barber is licensed under CC BY-SA 3.0. xのビット表現そのもの!
  10. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) 𝑥のビット表現の事を𝐵𝑥 と呼ぶことにする 𝐵𝑥 = 223𝐸 + 𝑀

    すると次のように近似で表現できる log2 (𝑥) ≈ 1 223 𝐵𝑥 + 𝜇 − 127 文章にすると、「F32の値xをintとしてスケー ルとバイアスをつけると𝐥𝐨𝐠𝟐 (𝒙) に近似でき る」となる。 "Float example" by Dave Barber is licensed under CC BY-SA 3.0.
  11. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) rsqrtの式をBを使って変形すると次のようになる 𝑦 = 1 √𝑥 log2 (𝑦)

    = log2 1 𝑥 1 223 𝐵𝑦 + 𝜇 − 127 ≈ − 1 2 log2 𝑥 1 223 𝐵𝑦 + 𝜇 − 127 ≈ − 1 2 1 223 𝐵𝑥 + 𝜇 − 127 𝑩𝒚 ≈ − 𝟏 𝟐 𝑩𝒙 + 𝟑 𝟐 𝟐𝟐𝟑(𝟏𝟐𝟕 − 𝝁) "Float example" by Dave Barber is licensed under CC BY-SA 3.0.
  12. 0x5f3759df-(bit_cast<uint32_t>(x) >> 1) 0x5f3759df - (bit_cast<uint32_t>(x) >> 1) 𝑩𝒚 ≈

    − 𝟏 𝟐 𝑩𝒙 + 𝟑 𝟐 𝟐𝟐𝟑(𝟏𝟐𝟕 − 𝝁)
  13. Magic number 青の線は次の式であった 𝑦 = log2 (𝑥 + 1) これをよく近似する次の式の𝜇を探索する

    𝑦 = 𝑥 + 𝜇 𝜇の取りうる範囲は次の通り [0, − log2 ln 2 − 1 ln 2 + 1] だいたい0.0から0.086くらいの値まで。 先ほどのBでいうところの、 0x5F2F796C から 0x5F400000
  14. Magic number 各Magic Numberに対して、 [10−10, 1010]の範囲で1万サンプルして相対誤差を出す x_values = np.logspace(-10, 10,

    num=100000, dtype=np.float32) 0x5F375A80 が最小になった。最大誤差: 1.74493e-03 [2]の0x5F375A86 は最大誤差: 1.75125e-03 (全floatの値をみたらこちらが勝つ)
  15. pow(x,?) 0x5f3759df - (bit_cast<uint32_t>(x) >> 1) 𝑩𝒚 ≈ − 𝟏

    𝟐 𝑩𝒙 + 𝟑 𝟐 𝟐𝟐𝟑(𝟏𝟐𝟕 − 𝝁) ↑「− 𝟏 𝟐 をかけてる」 = 「 − 𝟏 𝟐 乗してる」
  16. pow(x,?) -1/2乗以外に(おおよそ)任意の(ラフな)powが実装できる +(bit_cast<uint32_t>(x) >> 1) + M : ラフな平方根 +(bit_cast<uint32_t>(x)

    >> 2) + M : ラフな4乗根 -(bit_cast<uint32_t>(x) >> 2) + M : ラフな1/4乗根 +(bit_cast<uint32_t>(x)/5*11) + M : ラフなガンマ関数 +(bit_cast<uint32_t>(x) << 8) + M : ラフな128乗 かなり大雑把な近似で良く、極端な速度が求められるシチュエーションは何か?
  17. おわり float rsqrt(float x) { const uint32_t t = 0x5f3759df

    - (bit_cast<uint32_t>(x) >> 1); const float y = bit_cast<float>(t); return y * (1.5f - (x * 0.5f * y * y)); }
  18. 参考文献 [1]Wikipedia contributors. (2024, September 26). Fast inverse square root.

    In Wikipedia, The Free Encyclopedia. Retrieved from https://en.wikipedia.org/wiki/Fast_inverse_square_root [2]Leonid V. Moroz, Cezary J. Walczyk, Andriy Hrynchyshyn, Vijay Holimath, and Jan L. Cieslinski. 2016. Fast calculation of inverse square root with the use of magic constant - analytical approach. CoRR abs/1603.04483, (2016). Retrieved from http://arxiv.org/abs/1603.04483 [3]How the Quake engine cracked the inverse square root (2020, February 20). YouTube. Retrieved from https://www.youtube.com/watch?v=p8u_k2LIZyo