Making sense of floating point format + Q_sqrt from the Quake source

floating_point_sqrt_exploration.c

#include <stdio.h>

/*
  ruby -e '
     n = 0x5f3759df.to_s(2)
     n = "0 10000010 01011".ljust(32, "0").delete(" ")
     n = "0 01000001 00010000000000000000000"
     n = "0 10111110 01101110101100111011111"
     n = n.ljust(32, "0").delete(" ")
     sign     = (-1) ** n[0].to_i
     exponent = n[1,8].to_i(2)-127
     mantissa = 1 + n[9..].chars.each.with_index(1).map { |bit, exp| bit.to_f / (2**exp) }.sum
     result   = sign * (mantissa * (1 << exponent))
     formatted= "#{n[0]} #{n[1,8]} #{n[9..]}"
     pp n: n, sign: sign, exponent: exponent,
        mantissa: mantissa, result: result, formatted: formatted'
  {:n=>"01011111001101110101100111011111",
   :sign=>1,
   :exponent=>63,
   :mantissa=>1.4324301481246948,
   :result=>1.3211836172961055e+19,
   :formatted=>"0 10111110 01101110101100111011111"}
*/

void show_binary(char* desc, int32_t n) {
  int i = n;
  float f = * (float*) &n;
  char buffer[35];
  for(int i = 0; i < 35; ++i) buffer[i] = ' ';
  buffer[34] = '\0';

  // mantis, 23 bits
  for(int i = 0; i < 23; ++i) {
    buffer[33-i] = '0' + (n & 1);
    n >>= 1;
  }

  // exponent, 8 bits
  for(int i = 0; i < 8; ++i) {
    buffer[9-i] = '0' + (n & 1);
    n >>= 1;
  }

  // sign, 1 bit
  buffer[0] = '0'+n;

  printf("%s | %s (int: %d, float: %f)\n", desc, buffer, i, f);
}


// reciprocal sqrt
float Q_rsqrt(float number) {
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * (long*) &y;
  i  = 0x5f3759df - (i>>1);
  y  = * (float*) &i;
  y  = y * (threehalfs - (x2*y*y));
  return y;
}

int main() {
  /* printf("%d: %f\n", 9, Q_sqrt(9)); */

  const float number = 9;
  /* float number = 3; */

  /* 9  | 0 10000010 00100000000000000000000 */
  /* 3  | 0 10000000 10000000000000000000000 */

  int32_t i = * (long*) &number;
    show_binary("pre ", i);
    show_binary("shft", (i>>1));
    /* pre  | 0 10000010 00100000000000000000000 (int: 1091567616, float: 9.000000) */
    /* shft | 0 01000001 00010000000000000000000 (int: 545783808, float: 0.000000) */


  i = 0x5f3759df - (i>>1);
    show_binary("C   ", 0x5f3759df);
    show_binary("post", i);
    /* C    | 0 10111110 01101110101100111011111 (int: 1597463007, float: 13211836172961054720.000000) */
    /* post | 0 01111101 01011110101100111011111 (int: 1051679199, float: 0.342483) */

  float y = * (float*) &i; // y = C - (number/2)
  y = y * (1.5F - (number*0.5F*y*y));

  // (C-n/2) * (3/2 - n/2 * (C-n/2) * (C-n/2))
  // (C-n/2) * (3/2 - n/2 * (CC - Cn/2 - Cn/2 + nn/4))
  // (C-n/2) * (3/2 - n/2 * (CC - 2Cn/2 + nn/4))
  // (C-n/2) * (3/2 - n/2 * (CC - Cn + nn/4))
  // (C-n/2) * (3/2 - nCC/2 - nnC/2 + nnn/8)
  // (C-n/2) * (3 - nCC - nnC + nnn/4) * 1/2

  printf("%f\n", y);
  printf("%f\n", Q_rsqrt(number));
  /* show_binary( */

  /* return y; */
}