http://news.ycombinator.com/item?id=213056 <- Some comments

http://news.ycombinator.com/item?id=419166 <- Some comments

http://news.ycombinator.com/item?id=573912

http://news.ycombinator.com/item?id=896092 <- Some comments

http://news.ycombinator.com/item?id=1599635

http://news.ycombinator.com/item?id=2332793 <- Some comments

http://news.ycombinator.com/item?id=3115168 <- Some comments

http://news.ycombinator.com/item?id=3259199

This immediately suggests getting an initial guess for Newton's method by just shifting the bitwise form of an IEEE-754 floating point number by 1 and then doing some fix-ups for the negation.

The only hitch is that the exponent is stored in biased form rather than 2's complement form, so you have to first subtract out the bias, do the shift, and add back in the bias. Also, shifting can cause a least-significant 1 bit from the exponent to shift down into the most significant bit of the mantissa. But that turns out to be useful since it amounts to using the approximation sqrt(2) ~= 1.5. But doing all these bitwise operations is a lot of extra work compared to just doing a single shift. It turns out you can get almost the same effect with just a single addition of an immediate operand. Most of the bits of this addend are already predetermined by the need to handle the biasing and unbiasing in the exponent, and the remainder can easily be determined by brute force (the search space is only of size 2^24 or so) by just minimizing the approximation error over some test set. Doing that, I was able to converge to almost exactly the same 'magic number' as in Carmack's method.

My coworker Charles Bloom has a series of posts on this and related matters: http://cbloomrants.blogspot.com/2010/11/11-20-10-function-ap...

Here's an old email exchange I had with another coworker last time I was thinking about this stuff, where I tried to derive as many bits of the constant as possible by reasoning: https://gist.github.com/4e1d0dff97193fe5d745

Historical aside: How this algorithm came to be associated with Carmack is an amusing and roundabout tale. My partner in crime on Telemetry at RAD is Brian Hook. Hook was one of the first employees at 3dfx, where he designed and implemented the GLIDE graphics API. Gary Tarolli was the co-founder and CTO at 3dfx, and Hook got the rsqrt code from him. When Hook left 3dfx to work with Carmack on Quake 2 and 3 at id software, he brought that code with him and copied it into the codebase. Carmack never had anything to do with it. It's my recollection Tarolli didn't develop the algorithm either--he had gotten it from someone else in turn.

[0]:http://betterexplained.com/articles/understanding-quakes-fas...

> x = (float)&i;

breaks strict aliasing, and has undefined behavior. If you compile it with g++ -O2, you may get unexpected results. The solution is to use a union, or to compile with -no-strict-aliasing.

I used to have a lot of code like this, which worked fine on older compilers, but not on newer compilers. It's perfectly reasonable code, so I'm not sure why the newer compilers don't produce the expected behavior by default.

The reason is that there are optimizations that a compiler can do when it knows that two particular pointers will not reference the same location in memory at a given time (which is called aliasing). These optimizations have to do with memory access. If pointers cannot be proven not to alias each other, then the order in which they are loaded, modified, and stored will be very rigid since it is impossible to tell if changing the order will change the end result. Sometimes that's just life; you might have some complex interactions going with your pointers, and it may actually be that the order of operations can't be changed.

But if you can assume that the pointers don't alias each other, then the compiler is free to change the order of operations for better efficiency. It can often remove unnecessary loads and stores (if x and y don't alias, you don't have to reload y just because you changed x), and it can also put loads together at the beginning and stores together at the end (which is faster because, um, caching? I'm not really clear on that part).

The catch is that in general, there is no way for the compiler to know whether two pointers can alias each other. In C99, you can inform the compiler of this with the restrict keyword, but that requires that you do some careful analysis of your code to make sure you aren't lying to the compiler. If you mess up, you may end up with bugs that are very difficult to track.

So that's where strict aliasing comes in. Usually, two pointers of different types won't point at the same memory location. So they added the "strict aliasing" rule to make that fact official, because it means that in the vast majority of cases, code gets a nice speedup without any changes or difficult reasoning. Using a union is basically a way of informing the compiler that you need it to make an exception this one time.

I used to use a macro something like this to make the process concise:

  #define PUN(x, toType) (((union {__typeof__(x) __a; toType __b;})x).__b)

  int i = PUN(int, x);         // evil floating point bit level hacking
  i = 0x5f3759df - (i >> 1);  // what the fuck?
  x = PUN(float, i);
  x = x*(1.5f-(xhalf*x*x));

Modern CPU's can actually "look ahead" in the execution flow. So for something like this:

  MOV EAX,[memory_access]
  ADD EBX,ECX
  ADD EAX,EBX

Of course there are strange corner cases that have to be taken into account. For example, if the first instruction ends up causing an exception -- maybe the memory it refers to is in the swap file, or maybe it's dereferencing a NULL pointer -- the CPU has to be able to roll back the execution of any subsequent instructions, so it appears to the software as if the exception happened during the instruction that actually caused it.

If you're interested in the automatic program transformations performed by modern CPU's, you can read more by Googling for the Intel processor optimization manual [1] -- warning, PDF link. (I believe AMD publishes a similar document.)

[1] http://www.intel.com/content/dam/doc/manual/64-ia-32-archite...

  #define LOAD(TYP, ptr) (*((TYP *)memcpy((TYP[1]){ }, ptr, sizeof(TYP))))
  #define REINTERPRET(AS_TYP, FROM_TYP, val) LOAD(AS_TYP, (FROM_TYP[1]){ val })

> If the member used to read the contents of a union object is not the same as the member last used to store a value in the object, the appropriate part of the object representation of the value is reinterpreted as an object representation in the new type as described in 6.2.6 (a process sometimes called ‘‘type punning’’).

Please don't help spread the myth that compiler writers can break this idiom. That said, I use memcpy in my own code.

> When a value is stored in a member of an object of union type, the bytes of the object representation that do not correspond to that member but do correspond to other members take unspecified values.

So no undefined behaviour (nasal demons and such) but still unspecified.

union X{ char a; float b; }

If you store something in a, the other 3 bytes corresponding to b are unspecified

But I suspect this doesn't apply if a and b have the same size.

  import numpy
  import matplotlib.pyplot as plt
  import struct
  
  def power(x,p):
    i = struct.unpack("i", struct.pack("f", x))[0]
    i = (1-p) * 0x3f7a3bea + (p*i)
    return struct.unpack("f", struct.pack("i", i))[0]
  
  p = -0.5
  xs = numpy.arange(0.1,100,0.1)
  ys1 = map(lambda x: x**p, xs)
  ys2 = map(lambda x: power(x, p), xs)
  
  plt.plot(xs, ys1, label="Real function")
  plt.plot(xs, ys2, label="Approximation")
  
  plt.title("Graph of x^%d"%p)
  plt.legend(loc="upper left")
  plt.savefig("graph.png")

I'm guessing that (if I could find the article) the constant used in this floating point technique and the one used in this article are somehow inversely related too!

In any case, it's a very nice article ... thanks!

edit: http://blog.llvm.org/2011/05/what-every-c-programmer-should-...

Also because x * 1/sqrt(x) == sqrt(x) (so on a pentium, an sqrt was 3 cycles further an inv-sqrt, pretty cheap for an sqrt)