Square root union and Bit Shift

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?

The code uses a union to type-pun the float to an integer holding its bit-pattern. This is well-defined in C99 and later. (It's undefined behaviour in the C++ standard, but several mainstream C++ compilers also define the behaviour as an extension. This is a C question, but in C++ you'd use std::bit_cast<int32_t>(x)).

This doesn't depend on endianness, but does depend on int being the same width as float and float being IEEE-754 binary32, which is very common. Square root only works on non-negative numbers (sign bit clear), so int vs. unsigned doesn't matter, nor does 2's complement or arithmetic right shifts (implementation-defined behaviour of right-shifting a negative signed integer). Except for an input of -0.0, where this will shift the sign bit into the exponent field, producing a huge negative (or positive if using unsigned shifts) float.

---

To understand why this works, 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)

.