From mboxrd@z Thu Jan  1 00:00:00 1970
X-Spam-Checker-Version: SpamAssassin 3.4.5-pre1 (2020-06-20) on
	ip-172-31-74-118.ec2.internal
X-Spam-Level: 
X-Spam-Status: No, score=-0.5 required=3.0 tests=BAYES_05,FREEMAIL_FROM
	autolearn=ham autolearn_force=no version=3.4.5-pre1
X-Received: by 2002:ac8:6758:: with SMTP id n24mr67882915qtp.258.1609626390892;
        Sat, 02 Jan 2021 14:26:30 -0800 (PST)
X-Received: by 2002:a25:42c4:: with SMTP id p187mr94232630yba.504.1609626390636;
 Sat, 02 Jan 2021 14:26:30 -0800 (PST)
Path: eternal-september.org!reader02.eternal-september.org!news.gegeweb.eu!gegeweb.org!usenet-fr.net!proxad.net!feeder1-2.proxad.net!209.85.160.216.MISMATCH!news-out.google.com!nntp.google.com!postnews.google.com!google-groups.googlegroups.com!not-for-mail
Newsgroups: comp.lang.ada
Date: Sat, 2 Jan 2021 14:26:30 -0800 (PST)
Complaints-To: groups-abuse@google.com
Injection-Info: google-groups.googlegroups.com; posting-host=2601:193:4103:71a0:2d7b:8228:e7db:2ef2;
 posting-account=1tLBmgoAAAAfy5sC3GUezzrpVNronPA-
NNTP-Posting-Host: 2601:193:4103:71a0:2d7b:8228:e7db:2ef2
User-Agent: G2/1.0
MIME-Version: 1.0
Message-ID: <c7131827-fb1c-49e0-bb39-7bd9522aebdan@googlegroups.com>
Subject: Quick inverse square root
From: Matt Borchers <mattborchers@gmail.com>
Injection-Date: Sat, 02 Jan 2021 22:26:30 +0000
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
Xref: reader02.eternal-september.org comp.lang.ada:61013
List-Id: <comp.lang.ada>

I'm sure many of you have seen the Fast Inverse SquareRoot algorithm from t=
he open source Quake III engine.  I just encountered it a few days ago.  He=
re 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 =3D 1.5F;
   =20
       x2 =3D number * 0.5F;
       y =3D number;
       i =3D *(long *) &y;
       i =3D 0x5f3759df - ( i >> 1 );
       y =3D *(float *) &i;
       y =3D 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 :=3D a;
        i : UNSIGNED_32;
        for i'Address use y'Address;
    begin
        i :=3D 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 f=
rom 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 fo=
r Float case.

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

A second question I have is how to make this a generic for any Floating poi=
nt type.  I can only think that I have to provide three things: not only th=
e 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 :=3D a;
        i : U;
        for i'Address use y'Address;
    begin
        i :=3D 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_r=
ight".  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 pro=
grammer would require internal knowledge to make use of the generic.  Any t=
houghts on how to get the compiler to compute the magic number in the gener=
ic at compile time?

Thanks,
Matt