1

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;

}

JacekDuszenko
  • 136
  • 2
  • 11
  • That pointer cast does not magically convert `float` to `int`. It casts the pointer type and supposedly outputs that bit pattern in `int` format. The value of `int i` cannot be `100` (unless there is a rare matching example). Every line after `printf` is irrelavant. – Weather Vane Oct 14 '18 at 21:15
  • I obviously know that every line after printf is irrelevant. What does it mean that it outputs bit "pattern" in int format? How does such pattern look like? Where can I read more about bit patterns in particular types? – JacekDuszenko Oct 14 '18 at 21:17
  • So why post them? What you have is a 32-bit value which is interpreted differently in its context. – Weather Vane Oct 14 '18 at 21:18
  • 1
    This square-root kludge relies on the implementation using the IEEE-754 basic 32-bit binary floating-point format, one description of which is [here](https://stackoverflow.com/a/11963554/298225). The code reinterprets bytes in a way that is not supported by the current C standard. It could be corrected by using a union or by copying bytes. The particular value 0x5f3759df is the encoding for the floating-point value 13211836172961054720, which has an exponent of 63 and a significand of 0x1.6eb3be. – Eric Postpischil Oct 14 '18 at 21:20
  • 1
    Understanding the kludge requires an understanding of the floating-point format and certain mathematics. I am sure there are explanations on the web, possibly on Stack Overflow. On any hardware with a square-root instruction, the kludge is obsolete. – Eric Postpischil Oct 14 '18 at 21:22
  • @EricPostpischil Well I kind of already understood the mathematics behind the kludge but the thing that rankles me is this float to int and another way around pointer conversion – JacekDuszenko Oct 14 '18 at 21:24
  • [Here](https://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/) is one explanation. – Eric Postpischil Oct 14 '18 at 21:24
  • I mean, why did i get 1120403456 as the result of the first conversion and what this number means? – JacekDuszenko Oct 14 '18 at 21:26
  • Because it is not a value conversion but a pointer conversion. The `float` value is meaningless as `int` unless you want to examine the [constituent parts](https://en.wikipedia.org/wiki/Single-precision_floating-point_format). – Weather Vane Oct 14 '18 at 21:27
  • @WeatherVane So I got an integer sequence of numbers after converting the float pointer to integer pointer and dereferencing it. Where did those numbers come from? The question that would really help me understand it more would be: What do I have to do "by hand" to obtain 100 (my function argument) from this sequence 1120403456 – JacekDuszenko Oct 14 '18 at 21:30
  • The `int ` value in decimal makes no sense if you want to examine the `float` representation. Better in hex `42C80000` or even binary, please refer to the floating point format I linked in my previous comment, which you cannot have done within just 2 minutes. – Weather Vane Oct 14 '18 at 21:33
  • Replace `printf()` with code to output binary representation, see https://stackoverflow.com/questions/1024389/print-an-int-in-binary-representation-using-c – Severin Pappadeux Oct 14 '18 at 21:34
  • @WeatherVane: OP wants to examine the `float` representation. They have repeatedly told you they want to understand how 100 comes from 1120403456. The answer is that 1120403456 is the `float` representation for 100. – Eric Postpischil Oct 14 '18 at 21:38
  • @EricPostpischil my best advice to you is post an answer, since I have repeatedly tried to direct OP to the floating point format. – Weather Vane Oct 14 '18 at 21:44

3 Answers3

2

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

Kyrill
  • 2,963
  • 1
  • 8
  • 12
  • So If I convert 0x42c80000 to binary and split the resulting number according to https://en.wikipedia.org/wiki/Single-precision_floating-point_format (and by that I mean to sign bit, significand and exponent), will I get 100 after conversion to IEEE-754 ? – JacekDuszenko Oct 14 '18 at 21:46
  • Yes. You can use some online floating-point representation tool to see that the bit patterns are identical. `x` and `i` contain the same bit pattern. – Kyrill Oct 14 '18 at 21:56
  • Thank you for help – JacekDuszenko Oct 14 '18 at 22:40
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 2133−127 = 26. (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: + 26 • 1.5625, which equals 64•1.5625, which is 100.

If the exponent code were 0, it represents 2−126, 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.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
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.000000000000000000000002 (1.0 in decimal) to 1.111111111111111111111112 (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 = 2E - 127 × M

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

        √V = 2(E - 127)/2 × √M = 2E/2 - 127 × 263 × √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 × 263 × √2 × 1/√M

noting that 2127/2 = 263.5 = 263 × √2.

If we take the integer representation, and shift it one bit to the right, we effectively halve E. To multiply by 263, we need to add 63 to the exponent; 63×223. 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)2 (1597463007 in decimal, 0x5f3759df in hexadecimal, a single-precision floating-point approximation of √2127 ≈ 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.)

Nominal Animal
  • 38,216
  • 5
  • 59
  • 86