Q_sqrt

======================================================
    Calculating square root and inverse square root
    for 64-bit floating point numbers.
======================================================

1.  The basic constant values
  > L constant is the maximum value can be stored in
    the number's fraction ( for IEEE754 its 2^52 )
    * L = 4503599627370496
    
  > B constant is the maximum value can be stored in
    the number's exponent ( ITS SIGNED! ) ( for IEEE754 its (2^(11-1) - 1) )
    * B = 1023

2.  C_sqrt - magic number formula
  > For square root, magic number can be calculated according to the formula:
        ===========================================
        ==  magic = 0.5 * L * ( B - 0.0450465 )  ==
        ===========================================
    
  > So, the magic for C_sqrt will be:
    0|  uint64_t magic = (uint64_t)((0.5 * 4503599627370496 * (1023 - 0.0450465)))
    1|  printf("0x%llX\n", magic);  // 0x1FF7A3BEA91D9B00

3.  C_sqrt - apply formula
  > Now we need the formula to modify converted to int64_t double, so:
    0|  int64_t i = *((int64_t *)&num);
    1|  i = /* ---> */ (magic + (i >> 1)) /* <--- */;
    2|  num = *((double *)&i);

4.  C_sqrt - newton method
  > Just like that:
    0|  const double nhalf = (0.5 * num);
    .|  // .....
    ?|  num *= (0.5 + ((nhalf / num) / num));

5.  C_sqrt - implementation
    0|  double C_sqrt(double num)
    1|  {
    2|      const double nhalf = (num * 0.5);
    3|      int64_t i = *((int64_t *)&num);
    4|      i = (0x1FF7A3BEA91D9B00 + (i >> 1));
    5|      num = *((double *)&i);
    6|      num *= (0.5 + ((nhalf / num) / num));
    7|      num *= (0.5 + ((nhalf / num) / num));   // optional
    8|      num *= (0.5 + ((nhalf / num) / num));   // optional
    9|      return num;
   10|  }

6.  C_rsqrt - magic number
  > For square root, magic number can be calculated according to the formula:
        ===========================================
        ==  magic = 1.5 * L * ( B - 0.0450465 )  ==
        ===========================================
    
  > So, the magic for C_rsqrt will be:
    0|  uint64_t magic = (uint64_t)((1.5 * 4503599627370496 * (1023 - 0.0450465)))
    1|  printf("0x%llX\n", magic);  // 0x5FE6EB3BFB58D000

7.  C_rsqrt - apply formula
  > Now we need the formula to modify converted to int64_t double, so:
    0|  int64_t i = *((int64_t *)&num);
    1|  i = /* ---> */ (magic - (i >> 1)) /* <--- */;
    2|  num = *((double *)&i);

8.  C_rsqrt - newton method
  > Just like that:
    0|  const double nhalf = (0.5 * num);
    .|  // .....
    ?|  num *= (1.5 - (nhalf * (num * num)));

9.  C_rsqrt - implementation
    0|  double C_rsqrt(double num)
    1|  {
    2|      const double nhalf = (num * 0.5);
    3|      int64_t i = *((int64_t *)&num);
    4|      i = (0x5FE6EB3BFB58D000 - (i >> 1));
    5|      num = *((double *)&i);
    6|      num *= (1.5 - (nhalf * (num * num)));
    7|      num *= (1.5 - (nhalf * (num * num)));   // optional
    8|      num *= (1.5 - (nhalf * (num * num)));   // optional
    9|      return num;
   10|  }