Fast inverse square root algorithm के बारे में सब कुछ

(github.com/francisrstokes)

5 पॉइंट द्वारा GN⁺ 2024-06-03 | 2 टिप्पणियां | WhatsApp पर शेयर करें

Quake 3 से मशहूर हुआ fast inverse square root उस समय का एक performance solution था, जो 1 / sqrt(x) को float bits की reinterpretation और Newton-Raphson correction से तेजी से approximate करता था
इसका core यह है कि IEEE-754 32-bit float का integer bit pattern scaled और shifted log2(x) approximation की तरह handle किया जा सकता है
0x5f3759df - (i >> 1) दरअसल log2(x^-0.5) = -0.5 * log2(x) को integer shift और subtraction में बदलने वाला रूप है, और magic constant 3/2 * 2^23 * (127 - σ) से आता है
इसके बाद y = y * (1.5 - 0.5x * y * y) को एक बार चलाकर Newton-Raphson correction लगाया जाता है, और Quake code में दूसरा iteration comment out किया गया है
1999 में lighting और 3D vector normalization के लिए inverse square root की जरूरत प्रति सेकंड सैकड़ों से हजारों बार पड़ती थी, लेकिन modern hardware में dedicated floating-point processing की वजह से इसी trick की practical usefulness कम हो गई है

Quake code क्या करता है

float Q_rsqrt(float number) {
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = *(long*)&y;
  i  = 0x5f3759df - ( i >> 1 );
  y  = *(float*)&i;
  y  = y * ( threehalfs - ( x2 * y * y ) );

  return y;
}

यह function number के लिए inverse square root 1 / sqrt(number) का approximate value calculate करता है
सबसे famous हिस्सा है float value को long की तरह interpret करने के बाद 0x5f3759df - (i >> 1) bit manipulation करना
1999 में, जब Quake 3 रिलीज हुआ था, inverse square root धीमा और महंगा operation था, और lighting equations व normalization की जरूरत वाले 3D vector calculations में इसकी जरूरत प्रति सेकंड सैकड़ों से हजारों बार पड़ती थी
modern hardware में ऐसे calculations या तो CPU पर नहीं चलते, या CPU पर चलने पर भी बेहतर dedicated floating-point hardware की वजह से तेज होते हैं

IEEE-754 32-bit float representation

32-bit float तीन हिस्सों से बना होता है
- Sign: 1 bit, positive/negative बताता है
- Exponent: 8 bits, value किस range में है यह तय करता है
- Mantissa: 23 bits, उस range के अंदर position को linearly दिखाता है
सामान्य value को इस रूप में interpret किया जाता है

N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)

B = 127 biased exponent में इस्तेमाल होने वाला bias value है, और actual exponent e = E - B है
mantissa में सिर्फ m multiply नहीं होता, बल्कि 1 + m form इस्तेमाल होता है
- अगर m = 0 है तो 2^e
- अगर m = 1 के करीब हो जाए, तो अगले exponent range 2^(e+1) से ठीक पहले तक represent करता है
अगर exponent के सभी bits 0 हैं, तो यह sub-normal number है, और formula अलग हो जाता है

N = (-1)^S * 2^-126 * m

sub-normal, 0 और 0 के बेहद करीब छोटे numbers represent करने के लिए जरूरी होता है
अगर exponent के सभी bits 1 हैं, तो इसे special value के रूप में handle किया जाता है
- E = 255, M = 0 होने पर Infinity या -Infinity
- M != 0 होने पर NaN

float bits को integer की तरह देखने पर बनने वाला log relation

अगर float की internal representation को 32-bit integer की तरह देखें, तो इसे इस formula से दिखाया जा सकता है

I_x = 2^31 S + 2^23 E + M

inverse square root positive input के लिए होता है, इसलिए S = 0 रखने पर formula सरल हो जाता है

L = 2^23
I_x = L E + M

समान exponent range के अंदर mantissa position को linearly दिखाता है, लेकिन exponent बढ़ने पर mantissa steps की समान संख्या number line के ज्यादा बड़े interval को cover करती है
- E = 127, यानी e = 0, लगभग [1, 2) range है
- E = 128, यानी e = 1, लगभग [2, 4) range है
- दोनों ranges में mantissa steps की संख्या समान है, लेकिन दूसरी range दोगुनी चौड़ी है
इसी structure की वजह से float के raw bit pattern को integer की तरह देखने पर logarithmic relation दिखाई देता है

raw bits `log2(x)` का approximation हैं

float के bit pattern को integer I_x के रूप में interpret करें, तो इसे log2(x) के लिए piecewise linear approximation की तरह देखा जा सकता है
इस relation को नीचे दिए approximation से व्यक्त किया जाता है

log2(x) ≈ I_x / L - B

raw bit integer को mantissa size L = 2^23 से divide करके और exponent bias B = 127 घटाने पर log2(x) के करीब value मिलती है
mantissa interval के अंदर logarithm को linear approximation के रूप में handle किया जाता है

log2(1 + x) ≈ x + σ

σ approximation को adjust करने वाला tuning parameter है, और x [0, 1] range के अंदर exponent interval में position दिखाता है

inverse square root को log identity में बदलना

target नीचे दी गई value निकालना है

y = 1 / sqrt(x)

इसे exponential form में इस तरह बदला जाता है

y = x^-0.5

log identity apply करने पर inverse square root calculation नीचे दिए relation में बदल जाता है

log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)

float bits log2(x) approximation की तरह behave करते हैं, इसका इस्तेमाल करके x की integer bit representation I_x से y की integer bit representation I_y को सीधे approximate किया जा सकता है

I_y ≈ -0.5 I_x + 1.5 L (B - σ)

यही formula Quake code की core one-liner तक ले जाता है

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

i >> 1 integer bits को 1 bit right shift करता है, यानी 1/2 से multiply करने जैसा काम करता है
आगे वाला constant 0x5f3759df, 1.5 * L * (B - σ) के बराबर है

`0x5f3759df` constant की असलियत

अगर σ = 0 रखें, तो constant इस तरह calculate होता है

1.5 * 2^23 * 127 = 1598029824

इस value का hexadecimal representation 0x5f400000 है
यह Quake के actual constant 0x5f3759df से 566817 अलग है
इस difference से Quake code के corresponding σ value को calculate करें, तो यह मिलता है

σ = 377878 / 2^23
σ = 0.04504656

C code में वही constant इस तरह calculate किया जा सकता है

int32_t compute_magic(void) {
  double sigma = 0.0450465;
  double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
  int32_t i = expression;
  return i;
}

// -> 0x5f3759df

यहां double इस्तेमाल किया गया है, और integer conversion bit reinterpretation नहीं बल्कि normal casting है
यह σ value approximation optimize करने के लिए चुना गया था, लेकिन यह actual optimal value नहीं है और इसे किसने बनाया, यह भी पक्का नहीं है

यह सिर्फ simple hacking क्यों नहीं है

0x5f3759df - (i >> 1) float के raw bits को log approximation मानने के fact का इस्तेमाल करके inverse square root initial value बनाने वाला formula है
यह complex mathematical relation पर आधारित है, लेकिन execution step में सिर्फ shift और subtraction जैसे fast operations इस्तेमाल करता है
उस समय महंगे operations को प्रति सेकंड हजारों बार process करना पड़ता था, इसलिए यह approach hardware constraints के हिसाब से engineering design बन गई
हालांकि, यह algorithm सिर्फ normal float पर काम करता है
- sub-normal values में log2(1 + x) ≈ x + σ approximation assumption सही नहीं बैठता
- sub-normal में असल में 0 + x के करीब form आती है, इसलिए approximation टूट जाता है

Newton-Raphson correction से error कम करना

bit manipulation से मिली initial value काफी अच्छी होती है, लेकिन measurable error बचता है
नीचे वाली line approximation को काफी improve करती है

y = y * ( threehalfs - ( x2 * y * y ) );

यह line Newton-Raphson method apply करने वाला form है
inverse square root problem को Newton method के मुताबिक ढालने के लिए इसे नीचे दिए function की root खोजने की problem में बदलते हैं

f(y) = 1 / y^2 - x = 0

Newton method current approximation y_n से बेहतर approximation y_(n+1) इस तरह बनाता है

y_(n+1) = y_n - f(y_n) / f'(y_n)

f(y) = y^-2 - x का derivative यह है

f'(y) = -2y^-3 = -2 / y^3

division के बिना Newton correction formula

Newton formula को सीधे इस्तेमाल करें, तो कई floating-point divisions लगते हैं
इस algorithm के तेज होने की एक वजह floating-point division से बचना है
algebraically simplify करने पर यह division के बिना सिर्फ multiplication इस्तेमाल करने वाला form बन जाता है

y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)

Quake code में x2 = number * 0.5F से 0.5x पहले ही calculate किया जाता है, और फिर नीचे वाली line में इस्तेमाल होता है

y = y * ( threehalfs - ( x2 * y * y ) );

इस एक iteration के बाद maximum absolute error 0.175% है, और कई cases में error इससे कम होता है
original code में दूसरा Newton iteration है, लेकिन comment out किया गया है

// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

origin और related algorithms

यह algorithm John Carmack ने invent नहीं किया था, और exact origin 100% निश्चित नहीं है
related context के लिए Beyond3D का article linked है: The truth is the exact origin is not 100% certain
Chris Lomont ने logarithm approximation step में optimal sigma value खोजने की कोशिश पर paper लिखा: InvSqrt.pdf
CORDIC एक algorithm है जो floating-point के बिना सिर्फ addition और bit shifts से sine और cosine calculate करता है, और fast inverse square root से इसकी detailed method काफी अलग है
दोनों algorithms में समानता यह है कि वे mathematical observation को उस दौर के hardware constraints के हिसाब से efficiently apply करते हैं

2 टिप्पणियां

joyfui 2024-06-03

कभी-कभी फिर से सामने आ जाने वाला एक दिलचस्प कोड..ह

GN⁺ 2024-06-03

Hacker News टिप्पणियाँ

अगर कंप्यूटर 1999 के बाद बना है, तो आम तौर पर वह SSE instruction set सपोर्ट करता है, और इसमें _mm_rsqrt_ps है, जो एक बार में 4 inverse square roots को और तेज़ी से निकालता है: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
फिर भी, यहाँ बताई गई तकनीक अभी पूरी तरह बेकार नहीं हुई है। float/int conversion तेज़ है, लेकिन ऐसा hardware अब भी है जिसमें rsqrt, sqrt, pow, log instructions नहीं होते, और ऐसे operations को इस trick से approximate किया जा सकता है
- SSE का floating-point reciprocal instruction Intel और AMD पर थोड़ा अलग result दे सकता है, इसलिए अगर आप PCs के बीच deterministic result की उम्मीद करते हैं तो सिरदर्द हो सकता है: https://robert.ocallahan.org/2021/09/rr-trace-portability-di...
- दिलचस्प बात यह है कि SSE में सामान्य square root instruction भी है, लेकिन वह inverse square root से काफी धीमा है, इसलिए अगर precision में कमी स्वीकार कर सकते हैं तो sqrt(x) को x * 1/sqrt(x) के रूप में compute करना तेज़ पड़ता है
- सच तो यह है कि कंप्यूटरों का भारी बहुमत SSE तो दूर, i386/amd64 instruction set तक सपोर्ट नहीं करता, और emulation के बिना इसे सपोर्ट करने का अनुपात लगातार घट रहा है
  GPU instruction sets, ARM, RISC-V, AVR, PIC, 8051, FPGA आदि में approximate inverse square root operation अक्सर built-in होता है, लेकिन शायद उसे भी ऐसे ही algorithm से implement किया गया होगा
लेख पर थोड़ा nitpick करूँ तो, यह कहना सही नहीं है कि इस तरह की computation आजकल CPU पर नहीं होती। यह आम गलतफहमी है कि games या floating-point operations से भरे apps अपने सारे floating-point operations GPU को सौंपना चाहते हैं
असल में GPU को सौंपना सिर्फ बड़े और uniform workloads के लिए ही ठीक होता है। अगर one-off vector normalization कर रहे हैं, जैसे किसी object को दूसरे object की ओर देखने के लिए rotation matrix बनाना, तो उसे CPU पर ही रखना तेज़ है। GPU transfer time को छोड़ भी दें, तो single floating-point operation में CPU तेज़ होता है, क्योंकि GPU आम तौर पर कम clock पर चलता है और अपने ऊँचे FLOP count को parallelism से हासिल करता है
- यहाँ शायद GPU नहीं, FPU की बात की गई है। पुराने समय में FPU asynchronous calculation करता था, और आज उसे CPU में integrated हिस्सा माना जाता है
मैंने MMIX implementation लिखकर देखा, और यह assumption रखा कि original input value 2^-1021 से बड़ी है
अगर रुचि हो, तो Wikipedia पर भी इस function और उसके इतिहास की अच्छी explanation है: https://en.wikipedia.org/wiki/Fast_inverse_square_root
ऐसी कुछ चीज़ें इकट्ठी की हैं: https://github.com/ncruces/fastmath/blob/main/fast.go
इससे जुड़ा StackOverflow post भी है: https://stackoverflow.com/questions/32042673/optimized-low-a...
- मैं ऐसे techniques का collection शुरू करने के बारे में सोच ही रहा था, ताकि पहले बनाया हुआ late-80s style 3D engine फिर से लिख सकूँ; यह useful है
- fastmath package के benchmarks भी देखना चाहूँगा
nitpick करने का समय आ गया है। float formula में typo है, -1^S नहीं बल्कि (-1)^S होना चाहिए। पहला हमेशा -1 बन जाता है
raw bit pattern को interpret करना log का piecewise linear approximation है—यह explanation भी सटीक नहीं है। blue graph के data points के बीच की lines असल में मौजूद नहीं हैं, और कोई bit आधा 1 set नहीं हो सकता। बल्कि यह log के discrete version के ज्यादा करीब है, और जो data points सचमुच मौजूद हैं, यानी जहाँ red और blue lines मिलती हैं, वे literally scaled और shifted log जैसे ही हैं। बाकी लेख अच्छा है
- ठीक से समझ नहीं आया। बहुत छोटे 6-bit float को सोचें: sign 1 bit, exponent 2 bits, mantissa 3 bits हों तो [010000, 010111] interval में 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75 आते हैं
  लेकिन इन numbers के base-2 logs जो mantissa imply करते हैं, वे क्रमशः .0000000, .0010101, .0101001, .0111010, .1001010, .1011001, .1100111, .1110100 हैं, और पहले को छोड़कर float के 001, 010 आदि जैसे नहीं हैं। [2,4) interval के floats linear spacing में हैं, लेकिन corresponding log वैसा नहीं है, इसलिए लेख के मुताबिक float को log का piecewise linear approximation माना जा सकता है
- यह continuous piecewise linear approximation नहीं, बल्कि discrete piecewise linear approximation है। यह बात सही है कि blue line continuous नहीं है, लेकिन interpretation गलत है। blue graph सिर्फ कुछ intersection points नहीं, बल्कि x-axis पर बराबर दूरी पर रखे 256 individual points से बना है
  अगर पूरा graph हो, तो piecewise linear pattern में 2^32 choices होंगी, लेकिन original post ने जो बनाया है वह ऐसा full graph नहीं है। चूँकि लेख 32-bit integers और IEEE-754 32-bit float operations की बात करता है, इसलिए explanation में “discrete” छोड़ देना मुझे ठीक लगता है
बहुत से दिलचस्प concepts समझाने वाला अच्छा लेख है, लेकिन एक section में algebraic derivation हैरतअंगेज़ रूप से खराब है
“पहले form से इस form तक जाने के exact steps बहुत हैं, लेकिन completeness के लिए सब शामिल किए हैं” के बाद वाली derivation में अनावश्यक steps बहुत हैं, और कई sign errors हैं जो एक-दूसरे को cancel कर देते हैं। खासकर दूसरी line से तीसरी line पर जाते समय negative sign सही से distribute नहीं किया गया। अगर दूसरी line के बाद से शुरू करें, तो y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) से y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n) तक कहीं ज्यादा छोटा रास्ता है, और बीच के steps भी सही हैं। algebra समझने वाले व्यक्ति के लिए ये सभी obvious steps ही हैं
मशहूर code snippet का magic number optimal constant नहीं है। कोई दूसरा constant इस्तेमाल करने से relative error शायद लगभग 0.5% और कम हो सकता है
उस समय absolute optimum खोजना मुश्किल रहा होगा, लेकिन अब यह अपेक्षाकृत आसान है। मैं भी कभी इस rabbit hole में गया था और मेरे पास (1/x^2) और (1/x) के optimal magic numbers ढूँढने वाला Jupyter notebook है
- उस सवाल को explore करने वाले paper का link लेख के नीचे है
इस लेख में मेरे लिए सबसे दिलचस्प चीज़ “How Java's Floating-Point Hurts Everyone Everywhere” link था: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
लेखक William Kahan हैं, जिन्हें “Old Man of Floating-Point” के नाम से भी जाना जाता है: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
- विषय से अलग, लेकिन मैंने JAVAhurt PDF पढ़ना शुरू किया और typesetting भयानक है। ऐसा लगता है जैसे कोई TeX package इस्तेमाल हुआ है जो words के बीच की spacing को बहुत ज्यादा, वह भी असमान रूप से फैलाता है, और जैसे किसी दूसरे document को OCR करने के बाद extra spaces घुस गए हों
  fixed-width font वाले हिस्सों में भी अजीब extra spacing है। पढ़ने पर focus करना सचमुच मुश्किल था, और मुझे पता है कि वास्तव में ऐसा नहीं है, फिर भी यह लगभग किसी science crank manifesto जैसा लगा
पहले देखा हुआ यह video सचमुच बहुत अच्छा था: https://www.youtube.com/watch?v=p8u_k2LIZyo

Fast inverse square root algorithm के बारे में सब कुछ

Quake code क्या करता है

IEEE-754 32-bit float representation

float bits को integer की तरह देखने पर बनने वाला log relation

raw bits log2(x) का approximation हैं

inverse square root को log identity में बदलना

0x5f3759df constant की असलियत

यह सिर्फ simple hacking क्यों नहीं है

Newton-Raphson correction से error कम करना

division के बिना Newton correction formula

origin और related algorithms

संबंधित पढ़ाई

2 टिप्पणियां

Hacker News टिप्पणियाँ

raw bits `log2(x)` का approximation हैं

`0x5f3759df` constant की असलियत