@BJovke gives a nice answer that points at the truncated polynomial expansion as the key idea. This is a very standard approach in physics and maths. But still, why the magic number? The specific reason is the typical value of x. Let me explain.

The specific expansion being used is the https://en.wikipedia.org/wiki/Binomial_approximation . This is great for linearising small non-linear deviations from a default value. If you visit that Wiki page you'll find that you first need to write your real function in the form (1 + x)^a. But the inverse (really reciprocal) square root is x^(-1/2) so doesn't obviously have a suitable form for this -- there's no room for a + 1 offset on x.

Ah, but we could always write it as (b + x)^a then take out the b by division: b^a (1 + x/b)^a. You can also think of this as writing x in "natural units" of b.

For the reciprocal sqrt, a = -0.5, and so the first terms in the approximation are b^(-0.5) (1 - 0.5 x/b) = 1/sqrt(b) - 0.5 b^(-3/2) x. Note that this doesn't work for exactly b = 0, but 1/sqrt(x) is infinity for x = 0 so it's pretty reasonable that we can't actually go all the way to zero.

But, if you want to be valid in the vicinity of x = 0 then you want b as small as possible... i.e. b = 2^(-127) in a bit representation. Then that first term is sqrt(2^127) = 2^63.5. This makes sense: the default-best-guess term for reciprocal sqrt of a small x near that function's divergence at x = 0 is the sqrt of the largest representable number. And the long integer approximation to that, interpreted as a float, is the famous "magic constant". I've not checked, but I strongly suspect that the threehalfs term that follows it is implementing that b^(-3/2) in the second term of the approximation.

A corollary to this is that if you know that your reciprocal sqrts will be applied to numbers far from zero, you can probably get a better approximation by setting b to their typical value instead. I guess in Quake's raytracing, the x's tended to be small for whatever reason.