Re: Quick inverse square root

[Matt Borchers <mattborchers@gmail.com>]

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