Avoid Dereferencing in Carmacks Implementation of Approximate Roots

Just about everyone who has even looked into maths optimisation has come across the 'Quake III Magical Inverse Square Root Implementation' as used by John Carmack. For those that don't know it is a 32bit integer magic number implementation of the Newton–Raphson method for finding approximate roots.

The original code from Quake III looks like this:

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;
  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i>> 1 ); // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}

Okay, so what's the problem? The problem is the line commented as 'evil floating point bit level hacking' - the author didn't like it and I don't like it either. If you compile this code with gcc the following warning will be thrown:

warning: dereferencing type-punned pointer will break strict-aliasing rules

To find out what type-punning is and how to fix it we can turn to the gcc docs:

-fstrict-aliasing
Allows the compiler to assume the strictest aliasing rules applicable to the language being compiled. For C (and C++), this activates optimizations based on the type of expressions. In particular, an object of one type is assumed never to reside at the same address as an object of a different type, unless the types are almost the same. For example, an "unsigned int" can alias an "int", but not a "void*" or a "double". A character type may alias any other type.

Pay special attention to code like this:

            union a_union {
              int i;
              double d;
            };

            int f() {
              a_union t;
              t.d = 3.0;
              return t.i;
            }

The practice of reading from a different union member than the one most recently written to (called ``type-punning'') is common. Even with -fstrict-aliasing, type-punning is allowed, provided the memory is accessed through the union type. So, the code above will work as expected. However, this code might not:

            int f() {
              a_union t;
              int* ip;
              t.d = 3.0;
              ip = &t.i;
              return *ip;
            }

Every language that wishes to perform language-specific alias analysis should define a function that computes, given an "tree" node, an alias set for the node. Nodes in different alias sets are not allowed to alias. For an example, see the C front-end function "c_get_alias_set".

Enabled at levels -O2, -O3, -Os.

It is now clear what the problem is and writing my own implementation becomes trivial:

union _data32 {
  long l;
  float f;
} data32;
float Math::rsqrt( const float value ) {
    // Implementation of the Newton Approximation of roots.
    long i;
    float x, y;
    const float f = 1.5f;
    x = value * 0.5f;
    y = value;
    // i = *(long *) &y;
    data32.f = y;
    i = data32.l;
    i = 0x5f3759df - (i>> 1);
    // y = *(float *) &i;
    data32.l = i;
    y = data32.f;
    y = y * (f - (x * y * y));
    // Uncomment to calculate sqrt
    // y = y * (f - (x * y * y));
    return y;
}

- FIN -

210e

About this entry

You’re currently reading “Avoid Dereferencing in Carmacks Implementation of Approximate Roots,” an entry on ReloadSystems

Published:
11.08.09 / 12pm
Category:
Theory and Formulas, C++, 3D