Thanks for posting this, Sander.
There are other examples of bit hacks and I have always wondered whether others are out there waiting to be discovered. Here there is the possibility of long range effects due to carrying. More wild possibilities come from multiplication, of course counting on the conventions for overflow, as well as the integer encoding. The thinking from Wolfram's "A New Kind of Science" is relevant, and probably his methods as well.
Before I get too off topic, here are some definitions to explore this wholly within Wolfram Language:
magicconst = 16^^5f3759df;
FloatBitsToReal[list_] :=
N[(-1)^list[[1]] 2^(FromDigits[list[[2 ;; 9]], 2] - 127) FromDigits[
Prepend[list[[10 ;;]], 1], 2]/2^23];
RealToFloatBits[x_] := Join[{Sign[x] /. {1 -> 0, -1 -> 1}},
IntegerDigits[#[[2]] + 126, 2, 8], Rest[#[[1]]]] &[
RealDigits[x, 2, 24]];
and the C code then becomes
InvSqrt[x_] := Module[{i, x2, y},
x2 = .5 x;
y = x;
i = FromDigits[RealToFloatBits[y], 2];
i = magicconst - BitShiftRight[i, 1];
y = FloatBitsToReal[IntegerDigits[i, 2, 32]];
y*(1.5 - (x2*y*y))]
Here is an ArrayPlot where the first row is the magicconstant, the second row is the float for 1.2, then it is shifted to the right, and the last row is the subtraction from the magic constant.
I guess the reason why this works is that the base 2 exponent is effectively multiplied by -1/2 which serves as a pretty good approximation to x^(-1/2) no matter what happens to the mantissa. Simplifying the search for a magic constant to a manageable size, (akin to the NKS methodology) by looking for alternatives with the exponent part, we can just look at 8 bits.
MinimalBy[Range[0, 255],
Function[m,
Sum[Total[
IntegerDigits[
BitXor[127 - Quotient[i - 127, 2], m - BitShiftRight[i, 1]],
2]], {i, 0, 255}]]] (* {190, 191}*)
ArrayPlot[IntegerDigits[{190, 191}, 2, 9], Mesh -> True]
So (unless I made a mistake) there are two possibilities for the exponent part. The magic constant corresponds to 190. Maybe there is a better magic constant with the bits from 191?
