9 ms·
Fast Inverse Square Root
- timhh 6y agoThis is an explanation of the fast inverse square root that I hope is easy to understand. I also managed to improve on it a little bit.
- recursive 6y agoThanks for this. I have one nit. You show `q = 1598029824 - u/2;` as being identical to `q = 0x5F400000 - u >> 1;`, but every language I know of uses a different order of operations, giving different results. It might be clearest to provide parentheses in the second case.
- usmannk 6y agoInterestingly enough it seems to do the shift first in Swift
- TwoBit 6y agoAll well-known programming languages will do the shift first.
- thehobgoblin 6y agoI looked at about 20 different "popular" languages and I believe that only 2 (Go and Swift) of them have shifts at a higher precedence than addition. I think I'm missing a reference/joke =/
- hemdawgz 6y agoOr this might just be another confident but baseless HN commenter assertion.
- copascetic 6y agoC, C++, Java, and Rust all have shift as lower precedence than addition/subtraction. Edit: and Javascript, Python.
- recursive 6y agoYou can add C# as well.
- beached_whale 6y agoUse brackets to be explicit about your expected order of operations.
- timhh 6y agoOops! I added brackets, thanks!
- x3n0ph3n3 6y agoComputers do not have _real_ number systems, only _rational_ number systems and rational approximations of real numbers.
- jackpirate 6y agoThere are plenty of libraries for computable numbers [1], which sit between the rationals and reals [1] https://en.wikipedia.org/wiki/Computable_number https://en.wikipedia.org/wiki/Computable_number
- lifthrasiir 6y ago...which (naturally) can't determine if some number x is equal to y or not; that alone makes using constructive real numbers challenging and the usage is virtually non-existent. I believe the most-known software using them is the Android calculator since 6.0 and Hans Boehm (of the gc fame) documented various issues encountered [1]. [1] https://dl.acm.org/doi/fullHtml/10.1145/2911981 https://dl.acm.org/doi/fullHtml/10.1145/2911981
- longemen3000 6y agoNice, i was looking at this, but for sqrt (basically the same integer trick)
- User23 6y agoVery nice piece! Also, minuend and subtrahend are the standard terms.
- sika_grr 6y agoYou mention the wiki page is badly written, and I have the same feeling whenever I read about anything mathematics related on it. But... it's written by volunteers, so I'll take what I can get. Since you obviously have the talent for explaining things, do you think the wiki page could be edited and improved for clarity?
- timhh 6y agoI have thought about it, but I always worry about making large edits to Wikipedia pages that I'll put a lot of work into it and then someone will come along and just revert it. Maybe that fear is unfounded - I'll put it on my (very long) todo list!
- tomstuart 6y agoThis is a great post! The setup could be clearer: I found it a little difficult to follow (during “a real number x … this value which I will denote f … treat the 31-bit value as an unsigned integer u”) what each of the named variables was intended to refer to, and formatting the 31-bit number as 32 bits threw me off too. Nothing insurmountable but I could feel my brain stumbling before I got to the clever part.
- jedbrown 6y agoIt may be useful to mention that modern architectures often have vectorized instructions like vrsqrt14ps (accessible via the _mm512_rsqrt14_ps intrinsic) that provides a 14-bit approximation (there are more accurate variants too) in every lane with an inverse throughput of 2. These are faster than the integer bit hacks. https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=rsqrt14_ps&expand=4819 https://software.intel.com/sites/landingpage/IntrinsicsGuide...
- apocalypses 6y agoThough worth noting only for specific AVX-512 CPUs. Though you can still have vectorised support for regular precision for AVX2.
- jedbrown 6y ago_mm256_rsqrt_ps is 12-bit accuracy in AVX. https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=rsqrt_ps&expand=4819,4804 https://software.intel.com/sites/landingpage/IntrinsicsGuide...
- dheera 6y agoI imagine the bit hacks are still useful for microcontroller applications though.
- fluffything 6y agoI prefer architectures that have a vector instruction for computing one or two Newton iterations instead. That way you can quickly get the precision that you want ;)
- dragontamer 6y agovrsqrtps executes once-per-cycle throughput on Skylake / Icelake. You can always compute newton's method afterwards to improve the accuracy. But getting the maximum accuracy from one cycle is probably best.
- Greek0 6y agoHere's another article (from 2012) describing the same calculation: http://h14s.p5r.org/2012/09/0x5f3759df.html http://h14s.p5r.org/2012/09/0x5f3759df.html
- amatic 6y agoThe original article has the same error in the title. x^(-1/2) is the reciprocal square root. The fault is with math notation where rasing x to a negative power, i.e. x^-a means 1/(x^a); while inverse function F(x) is written F^-1(x). The inverse square root would be just square. Edit: I saw it called multiplicative inverse of the square root
- adamnemecek 6y agoGod, this is such an old meme at this point. What's next, telling me how good lisp is?
- IMTDb 6y agoObligatory XKCD: https://xkcd.com/1053/ https://xkcd.com/1053/
- adamnemecek 6y agoXKCD is legit shit. This one, the standards one, the Bobby Tables one. And the password one. Things that are barely worth a chuckle are repeated as gospel.
- deleted 6y ago[deleted]
- grayhatter 6y agoWow, bruh, who hurt you?
- nurettin 6y agoxkcd is as generic as they get. Randall goes idea -> generate numbers -> generate graph -> punchline and people milk it for as long as possible. It's like geek themed garfield. I wouldn't be surprised if he tries xkcd without stick people soon.
- grayhatter 6y agoSo it's generic and boring. But why are you really angry about it?
- adamnemecek 6y agoBecause it's everywhere.
- MauranKilom 6y agoSmall nitpicks: > But floating point numbers are always greater than LNS numbers Unless I'm missing something, they are never smaller, but they can be equal (at powers of two). > E.g. x1/2x^{1/2}x1/2, x−8x^{-8}x−8, x2x^2x2, though you probably wouldn't use it for positive exponents since you can just use multiplication to get an exact answer. 1/2 is positive. As is 1/3.
- timhh 6y ago> they can be equal Ha, funnily enough I did think of that, but I was trying to keep the explanation short and simple (and a bit hand-wavy). > 1/2 is positive. As is 1/3 Oops, fixed, thanks!
- echeese 6y agoI wonder how it holds up to something like this on a modern processor and compiler: float rsqrt(float number) { return 1.0f / sqrtf(number); }
- MauranKilom 6y agoWonder no longer! https://godbolt.org/z/M4TKb7 https://godbolt.org/z/M4TKb7 While this may look disappointing, it should be clear that adding 5% (or even 0.1%) error to inverse square root calculations is not something compilers are in the business of. Now, when you give the compiler more leeway (and -ffast-math is substantial leeway for anything less ephemeral than triangle normals on screen), you get much more interesting things: https://godbolt.org/z/85zG5r https://godbolt.org/z/85zG5r Turns out x86 has a (microcoded) instruction for this (and newer processors also have vectorized versions[1] of it). The compiler adds Newton-Rhapson for good measure. [1]: https://uops.info/html-lat/KBL/VRSQRTSS_XMM_XMM_XMM-Measurements.html https://uops.info/html-lat/KBL/VRSQRTSS_XMM_XMM_XMM-Measurem...
- SoSoRoCoCo 6y agoFYI: Type punning like this doesn't work on all compilers i = * ( long * ) &y; // evil floating point bit level hacking Learned this recently with Arm AC6. [Also, this kind of genius analytic approximation of hot functions that do math makes me all tingly]
- IAmLiterallyAB 6y agoWould a union do the trick?
- qppo 6y agoThat's the Linux way to do it union { int i; float f } u { .f = 1.23f }; int i = u.i; Another way is a memcpy, which I believe is the most defined way to do type punning int i; float f; memcpy(&i, &f, 4); But you also have to assume the size of those primitive types but that's pretty safe in modern C/C++.
- tsomctl 6y agoDo modern compilers optimize away the memcpy?
- qppo 6y agoyes
- colejohnson66 6y agoIn C++, type punning through pointers and unions is undefined behavior. Even `reinterpret_cast<T>` isn’t allowed because of aliasing (IIRC). The only “defined” way to do type punning is a memcpy. A compiler targeting something like x86 would optimize out the memcpy. For more information, see the C++20 final draft[0§7.6.1.9] [0]: https://isocpp.org/files/papers/N4860.pdf https://isocpp.org/files/papers/N4860.pdf
- 6y ago
- zeroonetwothree 6y agoTo me "inverse square root" means "square" so it makes the title of this kind of funny.
- drinkwell 6y agoYeah, shouldn't it be 'fast reciprocal square root'?
- MaxBarraclough 6y agoAgreed. Inverse is not a synonym for reciprocal, but it's too late, the name has stuck.
- jessriedel 6y agoNo, it's completely different than simple squaring. The inverse square root is undefined for negative numbers ;)
- MaxBarraclough 6y agoThat's not right. The so-called fast inverse square root function computes the reciprocal of the square root, not the square. You're not the only one to be confused by the name. An awful choice of name, in my opinion.
- balddenimhero 6y agoHere, "inverse" is an abbreviation of "multiplicative inverse" (aka. the reciprocal). Granted, shortening it to just "inverse" is misleading, but the name has stuck.
- person_of_color 6y agoWhat was the speedup from this? If only 1% of cpu time was spent on slow square root, and this sped it up by 100%, it would be barely worth it.
- taneq 6y agoI don't know how heavily it was used in Quake but I'd bet quite a lot. Carmack didn't mess around optimizing things for no reason, he was (and is) pretty focused on spending time (coder and CPU) where it will do most good. Also, a "mere" 1% speedup seems trivial for most general coding but squeezing 1% out of optimized game code is like blood from a stone. Add a bunch of similar tricks together and you can see 10-20% overall improvement which is huge.
- Kranar 6y agoI mean it's now ingrained in programming lore, but the reality is that Carmack has nothing to do with this and Carmack himself has never taken credit for it. The fast inverse square root dates back to the 80s and has its roots in SGI's systems. One of the developers who worked at SGI would eventually go on to work at ID developing Quake and worked on this optimization.
- taneq 6y agoOh, I didn't mean to imply Carmack invented it - he himself is pretty clear that he didn't. I just meant that incremental improvements in efficiency like this can add up surprisingly fast, and that once the low hanging fruit are taken, they're well worth doing.
- longemen3000 6y agoIt clearly depends on the workload. If you have tons of 3d articles or polygons on a screen, one of the most repeated operations are the calculation of norms (x/sqrt(sum(x)). John carmack came to the fast inverse sqrt root explicitly for this speed advantage
- gruez 6y ago
- rustybolt 6y agoA little self-promotion: I've written a blog post to answer the question that this one ends with (how to optimize divisions by constant integers): https://rubenvannieuwpoort.nl/posts/division-by-constant-unsigned-integers.html https://rubenvannieuwpoort.nl/posts/division-by-constant-uns... https://ridiculousfish.com/blog/posts/labor-of-division-episode-i.html https://ridiculousfish.com/blog/posts/labor-of-division-epis... and https://ridiculousfish.com/blog/posts/labor-of-division-episode-iii.html https://ridiculousfish.com/blog/posts/labor-of-division-epis... also explain this (and are probably better written than my blogpost).
- DoingIsLearning 6y agoOff-topic comment, I really like the clutter-free look of your site, with the latex to html generator. For anyone curious like I was: https://github.com/rubenvannieuwpoort/static-site-generator https://github.com/rubenvannieuwpoort/static-site-generator
- DarkWiiPlayer 6y agoI totally agree! It's easy to read and has a bit of a "paper" feel to it, but at the same time the sans-serif font makes it still feel "modern" and like a website.
- unwind 6y agoVery nice page! Minor nit, there's a typo here in the power-of-two example: uint divide(uint n) { return n << p; } the shifts should be to the right, obviously.
- rustybolt 6y agoThanks! I have a completely rewritten version in the pipeline, where I use easier versions of the theorems/proofs and fix some mistakes (in fact, the proof of lemma 1 is quite nonsensical if you look closely).
- ambar123 6y agoHow to find nth root of any real no...??
- colonwqbang 6y ago> [Fixed point arithmetic] has some niche uses but isn't commonly used It is very common in the embedded world and in hardware. I routinely program a CPU which has no FPU, so any arithmetic will have to be done in fixed point. It's perfectly possible to get work done with fixed point arith, it just requires a bit more thinking through so you don't run out of bits. Usually, I write unit tests where I compile my code for host architecture and compare the fixed point result to host's floating point. Then I can set a precise error margin for my XP approximation.
- DoingIsLearning 6y agoAnother 'niche' use case is doing DSP calculations on an FPGA. If you can accept the trade-off of some degree of error with fixed-point, you can have very fast and heavily parallelized implementations.
- nurettin 6y agoCopy of the original article that started all this https://www.beyond3d.com/content/articles/8/ https://www.beyond3d.com/content/articles/8/
- matsemann 6y ago> Games calculate square roots and inverse square roots all the time to find the lengths of vectors A trick here is that one often doesn't have to do the square root. For instance, if you want something to happen when an object is 5 units away from another object, it's normal to do if sqrt( (x2-x1)^2 - (y2-y1)^2 ) < 5 then ... but instead you can do if (x2-x1)^2 - (y2-y1)^2 < 5^2 then ... trading a sqrt for a squaring. And often the square of the distance can be cached or even a constant. So a lot of libraries (for instance vector library from libgdx[0]) contain a dst function, but also a dst2 function that skips the squaring to save some cycles when not needed. [0]: https://libgdx.badlogicgames.com/ci/nightlies/docs/api/com/badlogic/gdx/math/Vector2.html#dst2-com.badlogic.gdx.math.Vector2- https://libgdx.badlogicgames.com/ci/nightlies/docs/api/com/b...
- thdrdt 6y agoSometimes you can even get away by dropping the square completely: if (x2-x1) - (y2-y1) < (x4-x3) - (y4-y3) then ... I once used this in a path tracer to speed things up a little. The results where less accurate but sometimes this can be used as trade-off.
- TeMPOraL 6y agoI vaguely recall seeing in Graphics Gems an additional coefficient placed there that's been computed to minimize the error of this metric vs. Euclidean one.
- izackp 6y agohttps://gist.github.com/izackp/5ae96a740678946d8763c198b0e54d16 https://gist.github.com/izackp/5ae96a740678946d8763c198b0e54... I have a list of functions for distance calculation with metrics on speed and accuracy. If anyone would find it useful.
- DarkWiiPlayer 6y agoOne of my favourite hacks of all time (in part due to the comments)
- balddenimhero 6y agoIt is fun to read the source of old games for gems like the discussed `Q_rsqrt` function. However, one often wonders what the limits/assumptions/guarantees of such tricks are so I've blogged about SMT-based reasoning about such properties using `Q_rsqrt` as an example: https://bohlender.pro/blog/smt-based-optimisation-of-fast-inverse-square-root/ https://bohlender.pro/blog/smt-based-optimisation-of-fast-in... -- might be interesting to those not too familiar with the verification business.
- augustk 6y agoThe implementation can be more concisely written (and without reassigning variables) as float InvSqrt(float x) { long yl; float y; yl = 0x5f3759df - ((*(long *) &x) >> 1); y = *(float *) &yl; return y * (1.5F - (x * 0.5F * y * y)); }