Quick inverse square root

Quick inverse square root

[Matt Borchers]

I'm sure many of you have seen the Fast Inverse SquareRoot algorithm from the open source Quake III engine.  I just encountered it a few days ago.  Here it is, a bit reduced, from the original source:

    //C code from Quake III engine
    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;
    }

It is interesting how much clearer the Ada code version is:

    with Interfaces; use Interfaces;
    function QUICK_INVERSE_SQRT( a : FLOAT ) return FLOAT is
        y : FLOAT := a;
        i : UNSIGNED_32;
        for i'Address use y'Address;
    begin
        i := 16#5F3759DF# - shift_right( i, 1 );
        return y * (1.5 - (0.5 * a * y * y));
    end QUICK_INVERSE_SQRT;

The magic hexadecimal number is calculated from the formula:

  3/2 * 2**23 * (127 - mu) where mu is a constant close to 0.043.

My question is that I am trying to get this to work for Long_Float but I'm not having any luck.  I would expect that everything should be the same in the algorithm except for the types (Float -> Long_Float and Unsigned_32 -> Unsigned_64) and the "magic" hexadecimal number that should be calculated from the same formula but adjusted for the Long_Float bit layout.

    3/2 * 2**52 * (1023 - mu) where mu is the identical constant as used for Float case.

This doesn't seem to work and I haven't been able to find my error.  I'm sure it is something silly.  Does anybody have a suggestion?

A second question I have is how to make this a generic for any Floating point type.  I can only think that I have to provide three things: not only the obvious Float type, but also the Unsigned type of the same size, as well as the hex constant.

    generic
        type F is digits <>;
        type U is mod <>;
        magic : U;
    function G_Q_INV_SQRT( a : F ) return F;

I write the body like this:

    function G_Q_INV_SQRT( a : F ) return F is
        y : F := a;
        i : U;
        for i'Address use y'Address;
    begin
        i := magic - shift_right( i, 1 );
        return y * (1.5 - (0.5 * a * y * y));
    end G_Q_INV_SQRT;

    function QUICK_INVERSE_SQRT is
        new G_Q_INV_SQRT( FLOAT, UNSIGNED_32, 16#5F3759DF# );

This won't compile because the type U is not valid for the call to "shift_right".  How do I overcome this obstacle?

Once that is overcome, is there a way I can eliminate having to pass in the unsigned type along with the floating point type?  That seems like the programmer would require internal knowledge to make use of the generic.  Any thoughts on how to get the compiler to compute the magic number in the generic at compile time?

Thanks,
Matt

Re: Quick inverse square root

[Jeffrey R. Carter]

On 1/2/21 11:26 PM, Matt Borchers wrote:
> 
>      generic
>          type F is digits <>;
>          type U is mod <>;
>          magic : U;
>      function G_Q_INV_SQRT( a : F ) return F;
> 
> I write the body like this:
> 
>      function G_Q_INV_SQRT( a : F ) return F is
>          y : F := a;
>          i : U;
>          for i'Address use y'Address;
>      begin
>          i := magic - shift_right( i, 1 );
>          return y * (1.5 - (0.5 * a * y * y));
>      end G_Q_INV_SQRT;
> 
>      function QUICK_INVERSE_SQRT is
>          new G_Q_INV_SQRT( FLOAT, UNSIGNED_32, 16#5F3759DF# );
> 
> This won't compile because the type U is not valid for the call to "shift_right".  How do I overcome this obstacle?

Make it an explicit generic formal function parameter:

    with function Shift_Right (...) return ...;

> Once that is overcome, is there a way I can eliminate having to pass in the unsigned type along with the floating point type?  That seems like the programmer would require internal knowledge to make use of the generic.  Any thoughts on how to get the compiler to compute the magic number in the generic at compile time?

You would want to make use of the attributes of floating point types in ARM A.5.3

http://www.ada-auth.org/standards/rm12_w_tc1/html/RM-A-5-3.html

Whether these provide the information you need is another question. I don't see 
how you could declare the modular type in the generic.

-- 
Jeff Carter
"One day he told me he was a gynecologist. He
couldn't speak no foreign languages."
Take the Money and Run
147

Re: Quick inverse square root

[Dmitry A. Kazakov]

On 2021-01-02 23:26, Matt Borchers wrote:
> I'm sure many of you have seen the Fast Inverse SquareRoot algorithm from the open source Quake III engine.  I just encountered it a few days ago.  Here it is, a bit reduced, from the original source:
> 
>      //C code from Quake III engine
>      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;
>      }
> 
> It is interesting how much clearer the Ada code version is:
> 
>      with Interfaces; use Interfaces;
>      function QUICK_INVERSE_SQRT( a : FLOAT ) return FLOAT is
>          y : FLOAT := a;
>          i : UNSIGNED_32;
>          for i'Address use y'Address;
>      begin
>          i := 16#5F3759DF# - shift_right( i, 1 );
>          return y * (1.5 - (0.5 * a * y * y));
>      end QUICK_INVERSE_SQRT;

This is not equivalent to C code, you have likely a typo error.

> The magic hexadecimal number is calculated from the formula:
> 
>    3/2 * 2**23 * (127 - mu) where mu is a constant close to 0.043.
> 
> My question is that I am trying to get this to work for Long_Float but I'm not having any luck.  I would expect that everything should be the same in the algorithm except for the types (Float -> Long_Float and Unsigned_32 -> Unsigned_64) and the "magic" hexadecimal number that should be calculated from the same formula but adjusted for the Long_Float bit layout.
> 
>      3/2 * 2**52 * (1023 - mu) where mu is the identical constant as used for Float case.
> 
> This doesn't seem to work and I haven't been able to find my error.  I'm sure it is something silly.  Does anybody have a suggestion?

The formula you wrote above cannot be right. In effect, the factor y 
calculated from the exponent must be numerically same for both float 
(IEEE 754 single-precision floating-point) and double (IEEE 754 
single-precision floating-point). Which is apparently not. You should 
get the exponent multiplied by the same power of 2 as for float. For 
double, I make a wild guess, you should replace right shift by 1 with 
right shift by 30 = 32-2.

> A second question I have is how to make this a generic for any Floating point type.  I can only think that I have to provide three things: not only the obvious Float type, but also the Unsigned type of the same size, as well as the hex constant.
> 
>      generic
>          type F is digits <>;
>          type U is mod <>;
>          magic : U;
>      function G_Q_INV_SQRT( a : F ) return F;
> 
> I write the body like this:
> 
>      function G_Q_INV_SQRT( a : F ) return F is
>          y : F := a;
>          i : U;
>          for i'Address use y'Address;
>      begin
>          i := magic - shift_right( i, 1 );
>          return y * (1.5 - (0.5 * a * y * y));
>      end G_Q_INV_SQRT;
> 
>      function QUICK_INVERSE_SQRT is
>          new G_Q_INV_SQRT( FLOAT, UNSIGNED_32, 16#5F3759DF# );
> 
> This won't compile because the type U is not valid for the call to "shift_right".  How do I overcome this obstacle?
> 
> Once that is overcome, is there a way I can eliminate having to pass in the unsigned type along with the floating point type?  That seems like the programmer would require internal knowledge to make use of the generic.  Any thoughts on how to get the compiler to compute the magic number in the generic at compile time?

As Jeffrey already suggested, use S'Exponent to extract the exponent.

General notes.

1. C code relies on float being IEEE 754 single-precision floating-point 
number with endianness opposite to integer endianness numbers. The 
exponent must land in the integer's MSB. This is clearly non-portable.

2. The approximation is very crude. I am too lazy to estimate its 
precision within the intended range, which is what? [0, 1]?

3. Ergo, making it generic has no sense.

4. If you port it to Ada. Add assertions validating endianness and 
floating-point format.

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Quick inverse square root

[Matt Borchers]

On Sunday, January 3, 2021 at 5:58:41 AM UTC-5, Dmitry A. Kazakov wrote:
> On 2021-01-02 23:26, Matt Borchers wrote: 
> > I'm sure many of you have seen the Fast Inverse SquareRoot algorithm from the open source Quake III engine. I just encountered it a few days ago. Here it is, a bit reduced, from the original source: 
> > 
> > //C code from Quake III engine 
> > 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; 
> > } 
> > 
> > It is interesting how much clearer the Ada code version is: 
> > 
> > with Interfaces; use Interfaces; 
> > function QUICK_INVERSE_SQRT( a : FLOAT ) return FLOAT is 
> > y : FLOAT := a; 
> > i : UNSIGNED_32; 
> > for i'Address use y'Address; 
> > begin 
> > i := 16#5F3759DF# - shift_right( i, 1 ); 
> > return y * (1.5 - (0.5 * a * y * y)); 
> > end QUICK_INVERSE_SQRT;
> This is not equivalent to C code, you have likely a typo error.
> > The magic hexadecimal number is calculated from the formula: 
> > 
> > 3/2 * 2**23 * (127 - mu) where mu is a constant close to 0.043. 
> > 
> > My question is that I am trying to get this to work for Long_Float but I'm not having any luck. I would expect that everything should be the same in the algorithm except for the types (Float -> Long_Float and Unsigned_32 -> Unsigned_64) and the "magic" hexadecimal number that should be calculated from the same formula but adjusted for the Long_Float bit layout. 
> > 
> > 3/2 * 2**52 * (1023 - mu) where mu is the identical constant as used for Float case. 
> > 
> > This doesn't seem to work and I haven't been able to find my error. I'm sure it is something silly. Does anybody have a suggestion?
> The formula you wrote above cannot be right. In effect, the factor y 
> calculated from the exponent must be numerically same for both float 
> (IEEE 754 single-precision floating-point) and double (IEEE 754 
> single-precision floating-point). Which is apparently not. You should 
> get the exponent multiplied by the same power of 2 as for float. For 
> double, I make a wild guess, you should replace right shift by 1 with 
> right shift by 30 = 32-2.
> > A second question I have is how to make this a generic for any Floating point type. I can only think that I have to provide three things: not only the obvious Float type, but also the Unsigned type of the same size, as well as the hex constant. 
> > 
> > generic 
> > type F is digits <>; 
> > type U is mod <>; 
> > magic : U; 
> > function G_Q_INV_SQRT( a : F ) return F; 
> > 
> > I write the body like this: 
> > 
> > function G_Q_INV_SQRT( a : F ) return F is 
> > y : F := a; 
> > i : U; 
> > for i'Address use y'Address; 
> > begin 
> > i := magic - shift_right( i, 1 ); 
> > return y * (1.5 - (0.5 * a * y * y)); 
> > end G_Q_INV_SQRT; 
> > 
> > function QUICK_INVERSE_SQRT is 
> > new G_Q_INV_SQRT( FLOAT, UNSIGNED_32, 16#5F3759DF# ); 
> > 
> > This won't compile because the type U is not valid for the call to "shift_right". How do I overcome this obstacle? 
> > 
> > Once that is overcome, is there a way I can eliminate having to pass in the unsigned type along with the floating point type? That seems like the programmer would require internal knowledge to make use of the generic. Any thoughts on how to get the compiler to compute the magic number in the generic at compile time?
> As Jeffrey already suggested, use S'Exponent to extract the exponent. 
> 
> General notes. 
> 
> 1. C code relies on float being IEEE 754 single-precision floating-point 
> number with endianness opposite to integer endianness numbers. The 
> exponent must land in the integer's MSB. This is clearly non-portable. 
> 
> 2. The approximation is very crude. I am too lazy to estimate its 
> precision within the intended range, which is what? [0, 1]? 
> 
> 3. Ergo, making it generic has no sense. 
> 
> 4. If you port it to Ada. Add assertions validating endianness and 
> floating-point format. 
> 
> -- 
> Regards, 
> Dmitry A. Kazakov 
> http://www.dmitry-kazakov.de


Thank you Jeff and Dmitry.  I have a generic functioning now.

Jeff,
Using attributes I was able to come up with a magic number using:

magic : constant U := U(3.0 / 2.0 * 2.0**(F'Machine_Mantissa - 1) * (F(F'Machine_Emax - 1) - 0.043));

where U is the modular type of the generic and F the floating point type.  I find it quite unfortunate that the modular type needs to be exposed to instantiate the generic and hope to find a way to eliminate it.  I have not looked at the executable to see if the magic number has been pre-computed by the compiler and optimized away.

Also, I should have thought to pass 'shift_right' as a generic function parameter.  Thank you.

Dimitry,
I can't find the typo you are referring to nor how it is not equivalent and I think the math is correct.
I understand your points and its limited usefulness, but those limitations are also in the original C source and this was really only a curiosity project.

When people tell me that they use C for its low-level power and simplicity, like bit manipulations, and claim that other languages can't match C in that sense, I like to show them just how much better Ada can be -- aside from all the other benefits we all know.  Eliminating the generic, I think the main algorithm is much more clear in the Ada version.

Here's my final code which seems to work well enough on my machine.  The compiler required me to instantiate the generic with different names and then use renames for the function in the package specification. 

    with INTERFACES; use INTERFACES;

    generic
        type F is digits <>;
        type U is mod <>;
        with function SHIFT_RIGHT( n : U; amount : NATURAL ) return U;
    function G_QUICK_INVERSE_SQRT( a : F ) return F;

    function G_QUICK_INVERSE_SQRT( a : F ) return F is
        magic : constant U := U(1.5 * 2.0**(F'Machine_Mantissa - 1) * (F(F'Machine_Emax - 1) - 0.043));
        y : F := a;
        i : U;
        for i'Address use y'Address;
    begin
        i := magic - shift_right( i, 1 );
        return y * (1.5 - (0.5 * a * y * y));
    end G_QUICK_INVERSE_SQRT;

    function QINVSQRT is
        new G_QUICK_INVERSE_SQRT( LONG_FLOAT, UNSIGNED_64, shift_right );

    function QUICK_INVERSE_SQRT( a : LONG_FLOAT ) return LONG_FLOAT renames QINVSQRT;

    function QINVSQRT is
        new G_QUICK_INVERSE_SQRT( FLOAT, UNSIGNED_32, shift_right );

    function QUICK_INVERSE_SQRT( a : FLOAT ) return FLOAT renames QINVSQRT;

Re: Quick inverse square root

[Jeffrey R. Carter]

On 1/3/21 11:31 PM, Matt Borchers wrote:
> 
> Thank you Jeff and Dmitry.  I have a generic functioning now.

Glad to have been of help.

Regarding the unsigned type, it seems this only works if F'Size = 32 or 64, so 
you could write versions that use Unsigned_32 and Unsigned_64, and then make 
your generic function do

if F'Size = 32 then
    return QISR32 (A);
elsif F'Size = 64 then
    return QISR64 (A);
else
    raise Program_Error with "F'Size must be 32 or 64";
end if;

But I don't understand why this exists. In what way is it better than the 
(inverse) Sqrt operation of the FPU?

-- 
Jeff Carter
"I'm a kike, a yid, a heebie, a hook nose! I'm Kosher,
Mum! I'm a Red Sea pedestrian, and proud of it!"
Monty Python's Life of Brian
77

Re: Quick inverse square root

[Matt Borchers]

On Sunday, January 3, 2021 at 6:47:15 PM UTC-5, Jeffrey R. Carter wrote:
> On 1/3/21 11:31 PM, Matt Borchers wrote: 
> > 
> > Thank you Jeff and Dmitry. I have a generic functioning now.
> Glad to have been of help. 
> 
> Regarding the unsigned type, it seems this only works if F'Size = 32 or 64, so 
> you could write versions that use Unsigned_32 and Unsigned_64, and then make 
> your generic function do 
> 
> if F'Size = 32 then 
> return QISR32 (A); 
> elsif F'Size = 64 then 
> return QISR64 (A); 
> else 
> raise Program_Error with "F'Size must be 32 or 64"; 
> end if; 

This requires duplication of algorithm which is what the generic is supposed to avoid.  This would be like a  psuedo-generic.  I did try something similar for the magic number (putting both inside the generic with an if F'Size test), but the generic complained about the bit sizes being wrong in the opposite if branch for each of the instantiations and stated that a constraint error might be raised at run-time.  Of course the constraint error never occurred, but it was annoying to have the warnings and also annoying to have to disable the warnings with a pragma in the code to avoid them.

> But I don't understand why this exists. In what way is it better than the 
> (inverse) Sqrt operation of the FPU? 

I mentioned first that this code comes from the Quake III engine.  There must have been a purpose for it then or maybe it was never called but left in the source code.  There are many videos about it on YouTube.  I'm not really a low-level graphics guy, but I think it was intended to operate on the unit vector for intense graphics operations.

I think this algorithm would work on any floating point type with a bit layout similar to the IEEE-754 standard regardless of how many bits were allocated to the exponent and mantissa.

I don't have any personal use for it. It seemed like an easy example to show how Ada code can be simpler and just as powerful as C.  I tried to turn it into a generic just as an exercise in trying to eliminate the modular type from the generic interface after I realized that two types were required that were related only in bit size.  I did it with an if statement using F'Size (similar to what you did above), but it requires duplicating a lot of code because type declarations alone cannot be in if statements.  To avoid duplication, I think what would be necessary is a Type itself being a result of an expression.  Such as:

i : (if F'Size = 32 then Unsigned_32 else Unsigned_64);  

or a way to dynamically choose a modular type of a particular size based on input, such as:

i : Unsigned(F'Size)

Wait! Can I do it like this:

type Modular is mod 2**F'Size;

But then, of course, the built in shift functions won't work.  I need to investigate this option.

Another thought experiment now as I type this:  I think I can eliminate the modular type from the generic interface by declaring and overlaying a packed Boolean array on the generic float type and handling the bit shift myself.  Of course this runs counter to the original purpose of why this routine was created in the first place - speed of calculation.

Matt

Re: Quick inverse square root

[Matt Borchers]

I quickly found that F'Size is not static when F is a type used in a generic instantiation so the following does not work.

type Modular is mod 2**F'Size;

I then realized that the idea of using a packed Boolean array does not lend itself to easily transform the floating point magic number calculation to a packed Boolean array and then I would also have to deal with the subtraction as well as the shift.  So, this is also a non-starter.

Matt

Re: Quick inverse square root

[Jeffrey R. Carter]

On 1/4/21 4:50 AM, Matt Borchers wrote:
> On Sunday, January 3, 2021 at 6:47:15 PM UTC-5, Jeffrey R. Carter wrote:
>>
>> if F'Size = 32 then
>>    return QISR32 (A);
>> elsif F'Size = 64 then
>>    return QISR64 (A);
>> else
>>    raise Program_Error with "F'Size must be 32 or 64";
>> end if;
> 
> This requires duplication of algorithm which is what the generic is supposed to avoid.  This would be like a  psuedo-generic.

Well, no. Your original generic function would be inside the generic function, 
and QISR32 and QISR64 would be instantiations of it with Unsigned_32 and 
Unsigned_64 respectively.

-- 
Jeff Carter
"Why, the Mayflower was full of Fireflies, and a few
horseflies, too. The Fireflies were on the upper deck,
and the horseflies were on the Fireflies."
Duck Soup
95

Re: Quick inverse square root

[AdaMagica]

> with INTERFACES; use INTERFACES;
> generic 
> type F is digits <>; 
> type U is mod <>;
> with function SHIFT_RIGHT( n : U; amount : NATURAL ) return U; 
> function G_QUICK_INVERSE_SQRT( a : F ) return F; 
> 
> function G_QUICK_INVERSE_SQRT( a : F ) return F is 

  pragma Assert (F'Base'Size = U'Modulus);

> magic : constant U := U(1.5 * 2.0**(F'Machine_Mantissa - 1) * (F(F'Machine_Emax - 1) - 0.043));
> y : F := a; 
> i : U; 
> for i'Address use y'Address; 
> begin 
> i := magic - shift_right( i, 1 ); 
> return y * (1.5 - (0.5 * a * y * y));
> end G_QUICK_INVERSE_SQRT; 

Re: Quick inverse square root

[AdaMagica]

AdaMagica schrieb am Montag, 4. Januar 2021 um 12:13:53 UTC+1:
> > with INTERFACES; use INTERFACES; 
> > generic 
> > type F is digits <>; 
> > type U is mod <>; 
> > with function SHIFT_RIGHT( n : U; amount : NATURAL ) return U; 
> > function G_QUICK_INVERSE_SQRT( a : F ) return F; 

I haven't got the slightest idea for which range this function should be applied, but for sure not for the complete Float range.
Thus there should be an assertion about the range of the parameter a in the body. Or even better a precondition in the spec.

Re: Quick inverse square root

[Dmitry A. Kazakov]

On 2021-01-04 12:28, AdaMagica wrote:
> AdaMagica schrieb am Montag, 4. Januar 2021 um 12:13:53 UTC+1:
>>> with INTERFACES; use INTERFACES;
>>> generic
>>> type F is digits <>;
>>> type U is mod <>;
>>> with function SHIFT_RIGHT( n : U; amount : NATURAL ) return U;
>>> function G_QUICK_INVERSE_SQRT( a : F ) return F;
> 
> I haven't got the slightest idea for which range this function should be applied, but for sure not for the complete Float range.
> Thus there should be an assertion about the range of the parameter a in the body. Or even better a precondition in the spec.

It appears to be the Newton method with an heuristic used to choose the 
starting point. The descriptions is here:

    https://en.wikipedia.org/wiki/Fast_inverse_square_root

It also mentions a hack for double precision IEEE 754 floats.

P.S. The method makes no sense to implement or use on modern hardware.

-- 
Regards,
Dmitry A. Kazakov
http://www.dmitry-kazakov.de

Re: Quick inverse square root

[Egil H H]

For anyone interested, there's a discussion on the algorithm in this paper:

https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf

-- 
~egilhh

Re: Quick inverse square root

[Matt Borchers]

I never made any claims about the contemporary usefulness of this algorithm.  As I said, I turned a bunch of ugly C into a few lines of pretty Ada for a skeptic.  I'm sure there is no need for this now, but it is still historically interesting.  Trying to make a generic simply started as a thought experiment of how to eliminate the need for the two types (float and modular) which I now concede is impossible in its most generic form.

As computers get faster, storage gets larger, and code libraries get bigger, it is unfortunate that most programmers do not need to be as clever as they once were required to be.

Thanks for finding and sharing the PDF paper!  I'm amazed someone could write so many pages on this.  Out of curiosity, I need to compare my computed magic number with theirs from their 64-bit float example.

Matt

Re: Quick inverse square root

[Paul Rubin]

Matt Borchers <mattborchers@gmail.com> writes:
> As I said, I turned a bunch of ugly C into a few lines of pretty Ada
> for a skeptic.

You could write the C code similarly to the Ada, using a union.  I don't
know why the C program isn't written that way.  Maybe it is very old.

Re: Quick inverse square root

[Matt Borchers]

On Monday, January 4, 2021 at 4:06:03 PM UTC-5, Paul Rubin wrote:
> Matt Borchers <mattbo...@gmail.com> writes: 
> > As I said, I turned a bunch of ugly C into a few lines of pretty Ada 
> > for a skeptic.
> You could write the C code similarly to the Ada, using a union. I don't 
> know why the C program isn't written that way. Maybe it is very old.

That is a good point Paul.  At least Ada can "keep up" and do the low-level things that some C users claims to be best choice and thus the only good choice.

Re: Quick inverse square root

[Brian Drummond]

On Mon, 04 Jan 2021 12:55:33 -0800, Matt Borchers wrote:

> As computers get faster, storage gets larger, and code libraries get
> bigger, it is unfortunate that most programmers do not need to be as
> clever as they once were required to be.
> 
> Thanks for finding and sharing the PDF paper!  I'm amazed someone could
> write so many pages on this.  Out of curiosity, I need to compare my
> computed magic number with theirs from their 64-bit float example.

Having spent quite some time elsewhere getting sqrt down to a single 
clock cycle (throughput : 8 cycle latency) it doesn't surprise me at all. 
(The name Terje Mathisen comes to mind for assembly language 
implementations)

The odd coding (non use of union, strange use of intermediate variables) 
may well have been the result of compiler code generation limitations; 
the "better" form may have compiled to a few more instructions or run a 
little more slowly; not a good thing for a gamer on limited hardware!

Have you benchmarked the pretty Ada version against the original C ... or 
against a straightforward float operation on modern hardware?

-- Brian