Similar trick to that used in the fast inverse square root routine[1] popularized by Quake 3.
More detail on https://en.wikipedia.org/wiki/Fast_inverse_square_root#Alias... .
IEEE754 floats are a very well-designed binary format, and the fact that these approximations are possible is part of this design; indeed, the first known instance of this trick is for a=1/2 by W. Kahan, the main designer of IEE754.
If you want full precision, you need to do three Newton-Raphson iterations after the initial approximation. One iteration is:
y = y * (2.0F - x * y);