Origin of Quake3's Fast InvSqrt()

It's just a quick one page republish of something I'd written for my personal website back in the day, and linked here at some point or another, that we thought deserves some more limelight. One for the geeks (all of you then ) !

 

http://www.beyond3d.com/forum/showthread.php?t=28154

Can't get the staff...

Jawed

 

The article doesn't exist anymore on lungexpress (never rewrote my article display code), hence the republish and non-link to that thread, which I'm well aware exists

 

This article created a double flashback!

1) Buying a Deltron Realvision Flash3D Voodoo card in 1997 (my first 3D card).

2) Starting the mentioned NV40/Doom3 thread where DemoCoder posted the FastInvSqrt (first time I saw this code).

 

This 2003 paper offers a detailed analysis of the algorithm and its math, but doesn't seem to have much to offer on earlier history:

http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf

But as the response from Terje shows, the origin of the code is unlikely to have even been graphics-related. Chances are it wasn't even done by someone still programming. Serious floating point/fixed point assembly language hacking by uber math geeks has been going on since the fifties.

 

That paper is linked to in the article.

 

There are better ways

(double post)

 

It is faster and better to use the AltiVec/SSE instructions for reciprocal sqrt estimate. The reciprocal sqrt estimate instructions alone are more accurate than this method. Max error over the set of normal numbers:

this method: 28402 ulps
RSQRTSS: 4980 ulps
RSQRTSS + NR: 4 ulps (but Inv_sqrt(0) is NaN )

For altivec, you'd use vec_rsqrte instead. PowerPC also has a scalar reciprocal square root estimate instruction, which gets similar accuracy on some PowerPC machines (others get much less and need a NR iteration).

Note that none of these are a drop in replacement for division when stricter IEEE-754 conformance is required. Also note that in many cases, you don't need sqrt at all. For example, when looking to see if a point is inside a sphere, you can compare the square of the distance of the point from the center against the square of the radius.

Ian Ollmann

#include "xmmintrin.h"
//rsqrtss + Newton Raphson
float isqrt_est_nr( float x )
{
float xhalf = 0.5f * x;

//get reciprocal sqrt estimate
x = _mm_cvtss_f32( _mm_rsqrt_ss( _mm_set_ss( x ) ) );

x = x * (1.5f - xhalf*x*x);

return x;
}

//RSQRTSS
float isqrtf_est( float x )
{
return _mm_cvtss_f32( _mm_rsqrt_ss( _mm_set_ss( x ) ) );
}

 

Very true, but as the story says this code is ancient. There was no 3DNOW, MMX, SSE, AltiVec or anything back then. You either did an actual FP sqrt and divide (extremely slow), or you do an aproximation like this. No one would ever use something like this anymore, but it's interesting from a historical perspective.

Also good job geo slashdotting B3D . Glad you were only down for an hour or so.

 

B3D has not been down for me at all today, its been a little slow on loads but its been stable all the way.

 

Okay, trying 3rd time

(or are posts moderated before displaying?)

Did you ask <a href="http://trac.bookofhook.com/bookofhook/trac.cgi/wiki/AboutMe">Brian Hook</a>?

He seems like a likely candidate, or at least a conduit from 3dfx to id.

 

I'm not questioning you submitting the article, it's the type of thing I think /. should have more of. I was just commenting on the server's response. I read the story first then headed over to /. (yes, b3d is before /. on my morning reading list ) and saw it was posted. One of the tags was "slashdotted" so I tried coming back and the connection just hung.

 

Who wrote the code in HAKMEM?

The HAKMEM reference makes me think of Steve Russell, though I think I heard him say that Spacewar's numerics came from DEC.

However it's far more likely to be bill gosper; if not he would remember. Of course there was no IEEE format fp back then (there was barely an IEEE!) not to mention the 36-bit architecture of the PDP-6 with 18-bit floats so HAKMEM would really have just been a suggestion.

(don't forget the PDP-6 was the origin of BLT, the parent of bitblt, so no wisecracks about 36-bit words please. Those machines remain the most convenient and fun machines to program in assembly).

 

As it happens, the original author is now known ( ) after someone close to the code made themselves known, so we'll be running with a follow up piece tomorrow that explains all. Stay tuned!

 

02/23/2000 usenet post from Jason Dorie, in group comp.graphics.algorithms


http://groups.google.com/group/comp.graphics.algorithms/browse_frm/thread/a4bdef092236355a/3ec1b02e96dc6dbc?lnk=st&rnum=3#3ec1b02e96dc6dbc

float InvSqrt(float v)
{
float x2 = v * (float)0.5F;
float y = v;
long i = *(long *) &y;
i = 0x5f3759df - (i>>1);
y = *(float *)&i;

y = y * (1.5f - (x2 * y * y));
y = y * (1.5f - (x2 * y * y));
return y;

}



1997/10/13 post from Robert Harley

http://groups.google.com/group/comp.arch.arithmetic/browse_frm/thread/ddd30bfce90255b8/d2f08ddf951489fe?lnk=st&q=0x5fe80000UL&rnum=1#d2f08ddf951489fe

64 bit code, seems to use some of the same ideas.

 

Does anyone know if there is a version of this routine that doesn't use floats?