/* This is a novel and fast routine for the reciprocal square root of an
IEEE float (single precision). It was communicated to me by Mike Morton,
and has been analyzed substantially by Chris Lomont. Later (12/1/06)
Peter Capek pointed it out to me. See:

http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
http://playstation2-linux.com/download/p2lsd/fastrsqrt.pdf.
http://www.beyond3d.com/articles/fastinvsqrt/

The author of this has been researched but seems to be lost in history.
However, Gary Tarolli worked on it and helped to make it more widely
known, while he was at NVIDIA. Gary says it goes back to 1995 or
earlier. */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* ---------------------------- InvSqrt ----------------------------- */

float InvSqrt(float x) {
   float xhalf = 0.5f*x;
   int i = *(int *)&x;          // View x as an int.
   i = 0x5f3759df - (i >> 1);   // Initial guess.
   x = *(float *)&i;            // View initial guess as float.
   x = x*(1.5f - xhalf*x*x);    // Newton step.
   return x;
}

/* Notes: For more accuracy, repeat the Newton step (just duplicate the
line). The routine always gets the result too low. According to Chris
Lomont, the relative error is at most -0.00175228. Therefore, to cut its
relative error in half, making it approximately plus or minus 0.000876,
change the 1.5f in the Newton step to 1.500876f.
   Chris says that changing the hex constant to 0x5f375a86 reduces the
maximum relative error slightly, to 0.00175124. However, using that
value seems to usually give a slightly larger relative error.
   The routine can be adapted to IEEE double precision. */

/* ------------------------------ main ------------------------------ */

int main(int argc, char *argv[]) {
   float x = (float)atof(argv[1]);
   float yt = (float)(1.0/sqrt(x));
   float ya = InvSqrt(x);
   printf("x = %f\ntrue recip sqrt = %f, approx = %f, rel error = %f\n",
      x, yt, ya, (ya - yt)/yt);
   return 0;
}