Are we also shifting the mantissa in the fast inverse square root?

Question

I've read this blog about the algorithm:

http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

I understand that we have to manipulate the exponent in order to obtain the inverse square root. To this manipulation we take advantage from the IEEE 754 notation. My problem comes in this sentence (almost at the end) of the article:

"So, the code converts the floating-point number into an integer. It then shifts the bits by one, which means the exponent bits are divided by 2"

I we do that we are also manipulating the mantissa. Why we're doing this? As I said, the idea is to manipulate the exponent so, why we are also changing the mantissa? Can anybody help me? Thanks

Answer 1

You number is $x=2^{e-b}·(1+m)$ where $e$ is the exponent, $b$ the bias, $m$ the mantissa interpreted as fractional part. The inverse square root of this is $$ \frac1{\sqrt{x}}=2^{(b-e)/2}·(1-m/2+O(m^2)) $$ This shows you why the bit shift in the floating point number will come close to that result and also why you need to subtract it from a "magic" number, since the negated bit shift alone produces something like (but not exactly) $$ -2^{b+2-e/2}·(1-m/2) $$

---

To test this, compute the magic constant by reverse engineering the formula:

#include<stdio.h>
#include<math.h>

int main() {
    float x=1;
    int k;
    for(k=0; k<20; k++) {
        int i = *(int*)&x;
        float y = 1/sqrtf(x);
        int j = *(int*)&y;
        printf("x=%.20e : magic=%X\n",x,(i>>1)+j);
        x *= 1.06;
    }
    return 0;
}

with result

x=1.00000000000000000000e+00 : magic=5F400000
x=1.05999994277954101562e+00 : magic=5F3C7D3B
x=1.12359988689422607422e+00 : magic=5F396B7A
x=1.19101583957672119141e+00 : magic=5F36CCBA
x=1.26247680187225341797e+00 : magic=5F34A33A
x=1.33822536468505859375e+00 : magic=5F32F17E
x=1.41851890087127685547e+00 : magic=5F31BA4E
x=1.50363004207611083984e+00 : magic=5F3100C3
x=1.59384787082672119141e+00 : magic=5F30C841
x=1.68947875499725341797e+00 : magic=5F311480
x=1.79084753990173339844e+00 : magic=5F31E98F
x=1.89829838275909423828e+00 : magic=5F334BD6
x=2.01219630241394042969e+00 : magic=5F34DC35
x=2.13292813301086425781e+00 : magic=5F338AA4
x=2.26090383529663085938e+00 : magic=5F329A82
x=2.39655804634094238281e+00 : magic=5F320E46
x=2.54035162925720214844e+00 : magic=5F31E8A8
x=2.69277262687683105469e+00 : magic=5F322C9E
x=2.85433888435363769531e+00 : magic=5F32DD67
x=3.02559924125671386719e+00 : magic=5F33FE87

The magic constant used in the program, 5f3759df, is in the middle of that following some error minimization strategy or some other.

Answer 2

You need to know that in IEEE 754, the leading "1" in the mantissa is implicit -- it appears nowhere in the binary representation of the number. Consequently, if $M$ is the represented part of the mantissa and $e$ is the exponent, the represented number is $\left( 1+\frac{M}{2} \right) \times 2^{-e}$. As the other answer notes, $\sqrt{1+\frac{M}{2}} \approx 1+\frac{M}{4}$, which is the new mantissa when the represented part is shifted.