6 ms·
Everything I know about the fast inverse square root algorithm
- rogerallen 2y agoIf you're interested, Wikipedia also has a decent discussion of this function and it's history. https://en.wikipedia.org/wiki/Fast_inverse_square_root https://en.wikipedia.org/wiki/Fast_inverse_square_root
- schmorptron 2y agoI really liked this video about it when I saw it a while ago: https://www.youtube.com/watch?v=p8u_k2LIZyo https://www.youtube.com/watch?v=p8u_k2LIZyo
- jxyxfinite 2y agoDid John not write the code? Or did he copy paste it from somewhere?
- jb1991 2y agowikipedia says: > Brian Hook may have brought the algorithm from 3dfx to id Software.
- Rinzler89 2y agoAnd at 3dfx he probably learned it from some ex-SGI guy, and that guy learned it doing his PhD at Stanford from someone who worked for ed Ed Catmull and so on. The lore goes deep with stuff like this. There's rarely just one author but many who improve the formula over time.
- fragmede 2y agothere's an unrecognized genius out there that realized you can just apply a bit flip to get that out there. That John Carmack hasn't claimed it as his invention when he could have, speaks volumes about his character.
- Rinzler89 2y agoExcept he couldn't have even if he wanted to. There's always the risk some former 3dfx or SGI graybeard comes out with some paper binder with the original implementation, and calls you out on your bullshit if you do. When you're such a famous public figure no way you ever risk lying in public. Even if he were to getaway with it there was nothing for him to gain: he already has all the fame and money from honest work, there's no point risking it all to claim some piece of code.
- kqr 2y agoThen you call it a parallel discovery!
- Waterluvian 2y agoI’m not sure anyone should be impressed by someone not stealing even if they had a chance to. That’s just a basic human expectation.
- tsuica 2y agoYou'd think so.
- wongarsu 2y agoThe origin is really unclear. It's not certain which id employee added it to Quake 3 and where they got it from. John claims it wasn't him. Gary Tarolli, one of the founders of 3dfx, claims to remember tweaking the value of the hex constant to the value known today. But he says he didn't come up with the algorithm either. Maybe somebody else at 3dfx did, maybe he got it from somewhere else.
- pclmulqdq 2y agoI believe this particular line of code comes from Terje Mathisen, who is still active on the IEEE 754 standards committee today.
- johndough 2y agoIf your computer was built after 1999, it probably supports the SSE instruction set. It contains the _mm_rsqrt_ps instruction, which is faster and will give you four reciprocal square roots at once: https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#text=_mm_rsqrt_ps&ig_expand=5642&ssetechs=SSE https://www.intel.com/content/www/us/en/docs/intrinsics-guid... That being said, the techniques discussed here are not totally irrelevant (yet). There still exists some hardware with fast instructions for float/int conversion, but lacking rsqrt, sqrt, pow, log instructions, which can all be approximated with this nice trick.
- bnprks 2y agoAmusingly, (to me at least) there's also an SSE instruction for non-reciprocal square roots but it's so much slower than reciprocal square root that calculating sqrt(x) as x * 1/sqrt(x) is faster assuming you can tolerate the somewhat reduced precision.
- mabster 2y agoI wouldn't be surprised if _mm_rsqrt_ps is actually implemented using the same bit level trick. Same as Carmack's, we did a single step of Newton's method and it was definitely good enough.
- Findecanor 2y agoI dunno about Intel and AMD, but ARM and RISC-V use lookup tables for rsqrt. Unlike AMD and Intel, those tables are precisely defined in their respective specs.
- pclmulqdq 2y agoIntel provides bit-accurate code. In older chips it used a faithful bipartite ROM: https://web.archive.org/web/20120124193536id_/http://www.acsel-lab.com/arithmetic/papers/ARITH12/ARITH12_Das%20Sarma.pdf https://web.archive.org/web/20120124193536id_/http://www.acs...
- casey2 2y agoHe posts c code and says that it computes the inverse square root, but the code he posted is undefined c, so while it could compute the inverse square root of it's input it could also do a long list of other stuff. One of the many reasons you should stay away from real world languages when talking about algorithms unless you are an expert in the language.
- tadfisher 2y agoThe code he posted is straight from the Quake 3 source release, a program which has certainly been installed and used on hundreds on millions of real computers. UB pedantry makes no difference to the actual behavior of the program as demonstrated.
- wizzwizz4 2y agoIf we're talking real-world languages, the code is perfectly well-defined. It's undefined behaviour in Abstract C, the C of the Specification, which to my knowledge nobody has ever implemented, but Quake 3 was not written in Abstract C. It was written in Visual C++ 2003, gcc 2.95, and some other specific C implementation. The code snippet that the article claims is C: 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; } which, as far as I can tell, is perfectly well-defined Abstract C.
- deleted 2y ago[deleted]
- TillE 2y agoSpecifically, I think the original code will always do what you expect as long as sizeof(long) == sizeof(float), and the alignment is the same. In reality no compiler is gonna do weird stuff unless your hardware target is weird.
- kragen 2y agocompiler engineers have been disappointing our 'surely no compiler would ever be so perverse as to do x' expectations for 25 years now, and it doesn't seem that they're likely to stop soon
- atleastoptimal 2y agoSomething interesting is I've seen this mythologized as evidenced of the "cracked engineer" theory, wherein random engineers will accidentally stumble upon incredibly complex discoveries in the course of their day to day work and expect no fanfare for this, besting the scientists and researchers who take the traditional route but aren't spurred by necessity. People, like xkcd, purported this was either Carmack or some random employee who figured it out in a few hours when stuck on a problem then waltzed onto the next one. In truth, as noted in the document, it has a history that goes back decades in academia, this was simply the most notable time it was implemented. I think it speaks to an issue I see in the software engineering world where people assume that collaboration is for low-IQ people and all great innovation comes from some super-genius working on their own for long enough, and isn't required to share how they arrived at their findings. I'm sure the mythology of this algorithm is propelled somewhat by the enigmatic character of writing "what the fuck" next to adding the constant, implying a mystical element to its utility that was arrived at without needing to clarify it in some long boring research paper.
- galangalalgol 2y agoI think the xkcd comic was more about academia vs product engineering in for profit entities. If anything, it is the collaborative effort of many people working on a product with deadlines that ends up with some of them coming up with really interesting solutions just due to the number of people looking. The lack of sharing between companies does increase the incidence of the outsider effect, at the massive cost of reinventing basic knowledge over and over. The only reason breakthroughs still happen despite that is just the large number of people trying random stuff that probably won't work motivated by unrealistic deadlines. They don't know it probably won't work though...
- fragmede 2y agohttps://xkcd.com/664/ https://xkcd.com/664/ for the curious.
- atleastoptimal 2y ago
- layer8 2y agoWhat’s up with the diagonal lines in the graph backgrounds?
- zzo38computer 2y agoI wrote a implementation in MMIX: % Constants FISRCON GREG #5FE6EB50C7B537A9 THREHAF GREG #3FF8000000000000 % Save half of the original number OR $2,$0,0 INCH $2,#FFF0 % Bit level hacking SRU $1,$0,1 SUBU $0,FISRCON,$1 % First iteration FMUL $1,$2,$0 FMUL $1,$1,$0 FSUB $1,THREHAF,$1 FMUL $0,$0,$1 % Second iteration FMUL $1,$2,$0 FMUL $1,$1,$0 FSUB $1,THREHAF,$1 FMUL $0,$0,$1 This implementation makes an assumption that the original number is greater than 2^-1021.
- DrNosferatu 2y agoNote that you’re actually using 3 distinct magic numbers - not only the single 0x5f3759df, but also 0.5 and 1.5; So for further accuracy we can instead have: conv.i = 0x5F1FFFF9 - ( conv.i >> 1 ); conv.f = 0.703952253f ( 2.38924456f - x * conv.f * conv.f ); return conv.f; [from Wikipedia]
- ncruces 2y agoI collected a few of these here: https://github.com/ncruces/fastmath/blob/main/fast.go https://github.com/ncruces/fastmath/blob/main/fast.go Also see this StackOverflow: https://stackoverflow.com/questions/32042673/optimized-low-accuracy-approximation-to-rootnx-n https://stackoverflow.com/questions/32042673/optimized-low-a...
- qingcharles 2y agoThat's great, thank you. I was just thinking about starting a collection of these to begin rewriting my old 3D engine from the late 80s.
- nj5rq 2y agoPlease, post the progress if you end up rewriting it.
- koeng 2y agoI'd love to see some benchmarks on your fastmath package
- ncruces 2y agoIt's… probably not worth it. I just wrote this years ago (go.mod said 1.12) for the fun of it, thought I had it in a Gist/GitHub, and uploaded it yesterday in response to this post. One thing I remember trying was coding this up in ASM… which makes it worse, because it prevents inlining. But I learned the Go ASM syntax that way.
- nj5rq 2y agoVery interesting, thank you for posting it.
- qingcharles 2y agoI was building 3D engines a few years before the days of Quake, and having recently watched videos about optimizing the trig code in Super Mario 64, and these other faster algorithms, it makes me wonder how much more can be squeezed out of old hardware. I was optimizing my assembler to the nth degree, but optimizing the algorithm is always going to be the real winner.
- dahart 2y agoHere’s something new that isn’t mentioned and might be worth adding to the repo. (Edit I’m wrong, it is there, I just didn’t recognize it.) The magic number in the famous code snippet is not the optimal constant. You can do maybe 0.5% better relative error using a different constant. Maybe at the time it was infeasible to search for the absolutely optimal number, but now it’s relatively easy. I also went down this rabbit hole at some point, so I have a Jupyter notebook for finding the optimal magic number for this (1/x^2) and also for (1/x). Anyone wanna know what the optimal magic number is?
- bradleyjg 2y agoThere’s a link at the bottom to a paper exploring that question.
- dahart 2y agoGood point, I hadn’t read Lomont’s paper and should have. I read the section in the Wikipedia article talking about it, and did try the constant that it suggests, however it depends on doing extra Newton iterations and I looked at the relative error of the initial guess without Newton. I can see in the paper he found something within 1 bit of what I found. I’m not certain mine’s better, but my python script claims it is.
- qaisjp 2y ago> Anyone wanna know what the optimal magic number is? Sure.
- dahart 2y agoFor no Newton iterations, I thought I found 0x5f37641f had the lowest relative error, and I measure it at 3.4211. Of course I’m not super certain, Lomont’s effort is way more complete than mine. Lomont’s paper mentions 0x5f37642f with a relative error of 3.42128 and Wikipedia and the paper both talk about 0x5f375a86 being best when using Newton iterations.
- pcwalton 2y agoTo nitpick the article a bit: > It's important to note that this algorithm is very much of its time. Back when Quake 3 was released in 1999, computing an inverse square root was a slow, expensive process. The game had to compute hundreds or thousands of them per second in order to solve lighting equations, and other 3D vector calculations that rely on normalization. These days, on modern hardware, not only would a calculation like this not take place on the CPU, even if it did, it would be fast due to much more advanced dedicated floating point hardware. Calculations like this definitely take place on the CPU all the time. It's a common misconception that games and other FLOP-heavy apps want to offload all floating-point operations to the GPU. In fact, it really only makes sense to offload large uniform workloads to the GPU. If you're doing one-off vector normalization--say, as part of the rotation matrix construction needed to make one object face another--then you're going to want to stay on the CPU, because the CPU is faster at that. In fact, the CPU would remain faster at single floating point operations even if you didn't take the GPU transfer time into account--the GPU typically runs at a slower clock rate and relies on parallelism to achieve its high FLOP count.
- j16sdiz 2y agoI think he refers to FPU, not GPU. In the old days, FPU do async computation. FPU is now a considered an integrated part of CPU.
- solarized 2y agoNoob here. Please ELI 5 what this algo used in quake 3 ?. > Solve lighting equations Is it the light or laser effect when the gun is fired ?
- GrantMoyer 2y agoWhen you render a pixel on screen (of a diffuse object), to determine how brightly it's lit, you need to compute the angle between the normal at the point on the surface and the light source. Computing this angle requires normalizing one or more vectors, which is done by dividing by the magnitude, or the square root of the sum of the squares of the components. So the lights in this case are any lights in the scene, typically represented as a point of origin and a brightness.
- f4ncyp4ntz 2y ago[dead]
- f4ncyp4ntz 2y ago[dead]
- Exuma 2y agoI tried this in a raytracter today while learning rust and it made the scene render all white. Lol
- Jerrrry 2y agoThis actually made me chuckle, although it took an embarrassingly long second.
- g15jv2dp 2y agoTime for nitpicks, sorry. The formulas for the floats have typos. They should read (-1)^S, not -1^S (which always equals -1). Interpreting the raw bit patterns isn't a piecewise linear approximation of the logarithm. The lines between the data points on the blue graph don't actually exist, it's not possible for a bit to be half set to 1. It's rather a discrete version of the logarithm: the only data points that exist - where the red and blue lines meet - are literally equal to the (scaled, shifted) logarithm. Other than that, nice post!
- kqr 2y agoI'm not sure I understand. Imagine a very small 6-bit float, with 1 bit sign, 2 bits exponent, and 3 bits mantissa. The interval [010000, 010111] contains the numbers 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75. But! These numbers' base two logarithms imply mantissae of - .0000000 (equal to the float 000) - .0010101 (not equal to the float 001) - .0101001 (not equal to the float 010) - .0111010 (not equal to the float 011) - .1001010 (not equal to the float 100) - .1011001 (not equal to the float 101) - .1100111 (not equal to the float 110) - .1110100 (not equal to the float 111) because the floats in the [2,4) interval are linearly spaced, whereas the corresponding logarithms are not. In other words, the floats are a piecewise linear approximation of the logarithm – just as the article says.
- g15jv2dp 2y agoYou're right about the first part, of course. But don't you just need to truncate? In any case, my point was more that it's not "piecewise linear". Piecewise linear means that the map is defined on some interval; it's not, it's defined on a discrete set of values. To take your example, the map isn't defined on e.g., 2.3. You can decide to interpolate linearly between the value at 2.25 and the value at 2.5, but that's a decision you make which isn't reflected in the code. Or said differently: do you consider that any map defined on a discrete subset of R is really a piecewise linear map?
- kqr 2y agoIf your point is that it's more of an arithmetic sequence than a line segment, thanks to discretisation, then I agree but I don't have a good word for it other than "discrete piecewise linear". I mean it's not continuous, but it's still on an interval scale.
- p0seidon 2y agoThanks for stealing hours of my productive time going down this rabbit hole.
- timerol 2y agoThis is a good post explaining a lot of interesting concepts, with a section that has surprisingly bad algebra. > The exact steps to go from the first form to this one are numerous to say the least, but I've included it in full for completeness. The algebra included after that line has many unnecessary steps, as well as multiple cancelling sign errors. Most notably the second line to the third line does not distribute the negative sign correctly. Picking up after the second line, I would simply have: y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) y_n+1 = y_n + (1 - x * y_n^2) * (y_n / 2) y_n+1 = y_n + y_n / 2 - x * y_n^3/2 y_n+1 = 3/2 * y_n - x * y_n^3/2 y_n+1 = y_n (3/2 * y_n - x * y_n^2 / 2) y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n) This is fewer than 1/3 of the included steps (between the second line and the end), and has the advantage of being correct during the intermediate steps. I don't think any of my steps are anything other than self-evident to someone who understands algebra, but I'd be happy to entertain competing perspectives.
- omoikane 2y agoI thought the most interesting bit in this article was the link to "How Java's Floating-Point Hurts Everyone Everywhere": https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf Written by William Kahan, also known as the Old Man of Floating-Point: https://news.ycombinator.com/item?id=29042853 https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
- michaelcampbell 2y agoTotally not-related to the topic, but I started reading that JAVAhurt PDF and find the typesetting god-awful. Just me? Feels like they imported some TeX package to space out words on a line far, far too far... and inconsistently. Like it was OCR'd from something else and put in extra spaces. Even the monospace sections have weird extra spacing. Anyway, enough ranting for the day, but I found this really hard to focus on to read - almost felt like a science-kook type manifesto, though I know it is not.
- michaelcampbell 2y ago(post edit time edit...) I looked at some of his other writings and they all suffer this. I'm not sure what TeX template he's using, but it's bizarre. (IMO, of course.)