Square root union and Bit Shift

Question

I found this code with which the square root is obtained what surprises me is the way it does, using a union and bit shifts this is the code:

float sqrt3(const float x)  
{
  union
  {
    int i;
    float x;
  } u;

  u.x = x;
  u.i = (1<<29) + (u.i >> 1) - (1<<22); 
  return u.x;
} 

first is saved in u.x the value of x and then assign a value to u.i then the square root of the number and appears magically u.x

¿someone explain to me how this algorithm?

Answer 1

The above code exhibits UB (undefined behaviour), so it should not be trusted to work on any platform. This is because it writes to a member of a union and reads back from a member different from that it last used to write the union with. It also depends heavily on endianness (the ordering of the bytes within a multi-byte integer).

However, it generally will do what is expected, and to understand why it is worthwhile for you to read about the IEEE 754 binary32 floating-point format.

Crash Course in IEEE754 binary32 format

IEEE754 commonly divides a 32-bit float into 1 sign bit, 8 exponent bits and 23 mantissa bits, thus giving

                   31 30-23           22-0
Bit#:               ||------||---------------------|
Bit Representation: seeeeeeeemmmmmmmmmmmmmmmmmmmmmmm
Value:              sign * 1.mantissa * pow(2, exponent-127)

With the number essentially being in "scientific notation, base 2".

As a detail, the exponent is stored in a "biased" form (that is, it has a value 127 units too high). This is why we subtract 127 from the encoded exponent to get the "real" exponent.

Short Explanation

What your code does is it halves the exponent portion and damages the mantissa. This is done because the square root of a number has an exponent roughly half in magnitude.

Example in base 10

Assume we want the square root of 4000000 = 4*10^6.

4000000 ~ 4*10^6  <- Exponent is 6
4000    ~ 4*10^3  <- Divide exponent in half

Just by dividing the exponent 6 by 2, getting 3, and making it the new exponent we are already within the right order of magnitude, and much closer to the truth,

2000 = sqrt(4000000)

.