Re: Quick inverse square root

[Matt Borchers <mattborchers@gmail.com>]

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;