The voodoo has nothing to do with anything known to anyone in over 100 years ago.
The right shift not only does divide the mantissa in half but it adds 0.5 to it if the biased exponent is even it then multiplies the number by 2^-64
There was no mathematical representation that would allow doing all of that in the past.
Edit: plus there is an implied (1.) on the mantissa that is not represented on the shift, so it's (mantissa-1)/2+1
Then the subtraction goes further to subtract two numbers bit by bit with different biased exponents, so in theory what is doing is multiplying the number by 2^(numExp-ConstExp) doing the math, then shifting everything back and getting a good approximation that would cut on the iterations needed by using the newton method alone.
The constant is not arbitrary, I guess they started using log base 2 of 127, but then changed it to work with a specific range of numbers, seems the intention was to make it work with 55,300 in the middle but whoever did the program transposed the last two nibbles. But there is a number for which this approximation is tuned for, close enough to 55,300. I'll double check but I think it was 55,299.23 or something like that so it seems the intention wasn't targeted for that number. Edit: Correction, it waas 55299.32 I guess I have the same problem as Carmack
Yeah, if you plug 55,299.32 (0x47580352) that due to mantissa bits limitations is really 55299.32031, gives you the almost exact 1/sqrt(55299.32) = 0.00425246|12508714199066162109375, which is the closest representation in 32bit floats of the real answer (0.00425245870102170015628783736641)
The more you digress from that value (55299.32) the error increases, but by very little.
I'm surprised they didn't do the extra steps to use (n-55299.32) to get an extra compensation from the drift making it more accurate, and probably not needing the Newton method to bring the number more into convergence. Not a linear divergence (like 0x40 in the 65535 range and 0x20 in the 55300 range per integer increase) but seems like it would help get the solution closer with an extra multiplication so still no divisions.
Edit: What I meant by transposing the last two nibbles I mean the constant I think it was intended to be 5F3759FD which it will give the exact 1 over square root of 55300, maybe Carmack is a bit dyslectic.