What is the exact value of the first type-punning in fast square root inverse algorithm?

Question

Having this code, (this is a well-known fast square root inverse algorithm, that was used in Quake 3) I am unable to understand the print output. I understand the entire algorithm superficially but want to get an in-depth understanding. What is the value of i when printf prints it? This gives 1120403456 for me. Does it depend on the computer's architecture? I've read somewhere that type-punning like this results in undefined behaviour. On the other website I've read that the value of i at that point is the exact value of bits that are used by this variable. I am confused with this and honestly expected the value of i to be 100. How do I interpret the resulting 1120403456? How do I convert this value to decimal 100? Are those bits encoded somehow? This is the code excerpt:

#include<stdio.h>

int main()
{
float x = 100;
float xhalf = 0.5f*x;
int i = *(int*)&x;
printf("%d", i);
i = 0x5f3759df - (i>>1);
x = *(float*)&i;
x = x*(1.5f-xhalf*x*x);
return x;

}

Answer 1

It is much older than Quake III, although the exact magic constant has varied a bit.

The encoding of a finite single-precision floating point number according to IEEE-754 (which is what most modern computers use) is

     31 30           23 22                                          0
     ├─┼─┬─┬─┬─┬─┬─┬─┬─┼─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┤
     │S│   Exponent    │                   Mantissa                  │
     └─┴───────────────┴─────────────────────────────────────────────┘

The highest bit, S is the sign. The value is negative if it is set, positive otherwise.

If Exponent and Mantissa are both zero, then the value is either 0.0 (if S is zero) or -0.0 (if S is 1).

If Exponent is zero, but Mantissa nonzero, you have a denormal number very close to zero.

Representations where Exponent is 255 (all bits set) are reserved for infinity (if Mantissa is zero), and not-a-number (if Mantissa is nonzero).

For Exponent values from 1 to 254, inclusive, Mantissa contains the fractional bits of the mantissa, with an implicit 1 in the integers position (This is what the (1) means in e.g. 0 10000001 (1)1011011011011011011011.) In other words, the value Mantissa represents for these Exponent values, is from 1.00000000000000000000000₂ (1.0 in decimal) to 1.11111111111111111111111₂ (1.99999994 in decimal), inclusive.

Now, if we consider only nonnegative values (with S clear), with E = Exponent (between 1 and 254, inclusive), and M the logical value of Mantissa (between 1.0 and 1.99999994), then the value V the single-precision floating point represents is

        V = 2^{E - 127} × M

The 127 is the exponent bias. The square root of V is

        √V = 2^{(E - 127)/2} × √M = 2^{E/2 - 127} × 2⁶³ × √2 × √M

and the inverse square root, which is the target operation here,

        1/√V = 2^{(127 - E)/2} × 1/√M = 2^{-E/2 - 127} × 2⁶³ × √2 × 1/√M

noting that 2^127/2 = 2^63.5 = 2⁶³ × √2.

If we take the integer representation, and shift it one bit to the right, we effectively halve E. To multiply by 2⁶³, we need to add 63 to the exponent; 63×2²³. However, for the square root operation, instead of multiplying by √2 × √M, we effectively just multiply by M/2. This means that that (shifting the integer representation one bit right, then adding 63 to the exponent) yields an estimate of a square root that is off by a factor between 0.7071067 and 0.75 for arguments greater than 1.0, and by a factor between 0.7071067 and 1448.155 for values between 0 and 1. (It yields 5.421×10^-20 for √0, and 0.75 for √1.)

Note that for the inverse square root operation, we wish to negate E.

It turns out that if you shift the integer representation right by one bit, and then subtract that from (1 01111100 1101110101100111011111)₂ (1597463007 in decimal, 0x5f3759df in hexadecimal, a single-precision floating-point approximation of √2¹²⁷ ≈ 13043817825332782212), you get a pretty good approximation of the inverse square root (1/√V). For all finite normal arguments, the approximation is off by a factor of 0.965624 to 1.0339603 from the correct value; a very good approximation! All you need is one or two iterations of Newton's method for the inverse square root operation on top. (For f(X) = 0 iff X = 1/√V you need f(X) = 1/X^2 - V. So, each iteration in X is X - f(X)/f'(X) = X × (1.5 - 0.5 × V × X × X. Which should look quite familiar.)

Answer 2

To interpret 1120403456 or any other number as an IEEE basic 32-bit binary floating-point number:

• Write the number as 32 bits, in this case, 01000010110010000000000000000000.
• The first bit is the sign. 0 is + and 1 is −.
• The next eight bits are an encoding of the exponent. They are 10000101, which is 133 in decimal. The exponent is encoded with a bias, 127, meaning the actual power of two represented is 2^133−127 = 2⁶. (If the exponent is 0 or 255, it is a special value that has a different meaning.)
• The remaining 23 bits encode the significand. For normal exponents (codes 1 to 254), they encode a significand formed by starting with “1.” and appending the 23 bits, “10010000000000000000000”, to make a binary numeral, “1.10010000000000000000000”. In decimal, this is 1.5625.
• The full value encoded is the sign with the power of two multiplied by the significand: + 2⁶ • 1.5625, which equals 64•1.5625, which is 100.

If the exponent code were 0, it represents 2⁻¹²⁶, and the significand would be formed by starting with “0.” instead of “1.”

If the exponent code were 255, and the significand bits were zero, the object would represent infinity (+∞ or −∞ according to the sign bit). If the exponent code were 255, and the significand bits were not zero, the object would represent a Not a Number (NaN). There are additional semantics to NaN not covered in this answer.

Answer 3

The value of i printed after int i = *(int*)&x; is the bit representation of the floating-point 100.0, which is what x was initialised to. Since you're using the %d format in printf it will print that representation as a decimal integer.

The bitpattern of 100.0 in an IEEE 32-bit float is 0x42c80000, which is 1120403456 in decimal