3 ms·
I would probably write it this way today (well, actually I wouldn’t write it today, I would instead use the reciprocal square-root instructions provided by all
by stephencanon 14y ago
I would probably write it this way today (well, actually I wouldn’t write it today, I would instead use the reciprocal square-root instructions provided by all the common architectures, but let’s pretend):
static uint32_t toRepresentation(float x) {
uint32_t rep;
memcpy(&rep, &x, sizeof x);
return rep;
}
static float fromRepresentation(uint32_t rep) {
float x;
memcpy(&x, &rep, sizeof x);
return x;
}
static float rsqrt_linearApproximation(float x) {
uint32_t xrep = toRepresentation(x);
uint32_t yrep = 0x5f3759df - xrep/2;
return fromRepresentation(yrep);
}
static float rsqrt_newtonRaphsonStep(float x, float y) {
return y*(1.5f - (0.5f*x)*y*y);
}
float Q_rsqrt(float x) {
float y = rsqrt_linearApproximation(x);
y = rsqrt_newtonRaphsonStep(x, y);
y = rsqrt_newtonRaphsonStep(x, y);
return y;
}
The only semi-cryptic thing about it, to someone versed in numerics, is line with the magic number 0x5f3759df. While I would expect a professional to be able to quickly re-derive it and understand the intention, I would still write a comment, along the lines of “We approximate the floating-point representation of 1/sqrt(x) with a first-order minimax polynomial in the representation of x.”
All that said, I tend to write a lot of comments. Much of the code I write professionally is math library code, and it’s not uncommon for me to have a ratio of several paragraphs of error analysis or correctness proofs to a few lines of code.
- InclinedPlane 14y agoIt's not just semi-cryptic, it's magic. Your example above includes "xrep/2" instead of i>>1, but this does nothing to increase the understanding of what the fuck is going on here. This is a floating point number and some deep juju bit twiddling is happening with it, it's not just being divided by 2. The exponent is being effectively divided by 2 and the fractional part of the mantissa is also being effectively divided by 2. Except in the case where the representation of the exponent ended in a 1, in which case that 1 is shifted down into the fractional part of the mantissa. Exactly why this is ok and doesn't hurt the approximation is a subject that could fill a report.
- stephencanon 14y agoFloating point representations are approximately logarithms. It should be clear to anyone with working knowledge of floating point numerics that dividing by two and negating the representation of a floating point number corresponds to taking the reciprocal square root of the number. That's not magic, that's math. The constant simply accounts for the bias in the floating point representation and balances the error over the domain of interest. Again, not magic. There is no "bit-twiddling" at all.
- lmkg 14y agoExplaining the details of the bit-twiddles does not make it any less bit-twiddling. And explaining the details behind the opaque steps does not make them any less opaque to the uninitiated. You're using bit operations on a data type with a non-trivial bit structure. That's pretty much the definition of bit-twiddling. The fact that blunt bit operations on such a finely-defined structure can correspond to useful mathematical operations is non-obvious.
- mrbrowning 14y agoThe fast inverse square root trick is clever, but this is right on. It's bypassing the provided interface for FP operations and mucking around with implementation details. That's not to say that that's a bad thing, or that floating point numbers aren't a leaky abstraction, or that the real world doesn't demand this sort of optimization from time to time. But it's still true.
- stephencanon 14y ago> You're using bit operations on a data type with a non-trivial bit structure. It’s using integer operations on an encoding that has a meaningful (approximate) interpretation as an integer.
- coldtea 14y ago>And explaining the details behind the opaque steps does not make them any less opaque to the uninitiated. Emm, by definition it does.