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=-1.9 required=3.0 tests=BAYES_00,FREEMAIL_FROM autolearn=ham autolearn_force=no version=3.4.5-pre1 X-Received: by 2002:a37:5dc3:: with SMTP id r186mr28457182qkb.258.1609713076120; Sun, 03 Jan 2021 14:31:16 -0800 (PST) X-Received: by 2002:a5b:107:: with SMTP id 7mr106640754ybx.253.1609713075864; Sun, 03 Jan 2021 14:31:15 -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: Sun, 3 Jan 2021 14:31:15 -0800 (PST) In-Reply-To: Complaints-To: groups-abuse@google.com Injection-Info: google-groups.googlegroups.com; posting-host=2601:193:4103:71a0:e481:b83e:1f1a:a8ca; posting-account=1tLBmgoAAAAfy5sC3GUezzrpVNronPA- NNTP-Posting-Host: 2601:193:4103:71a0:e481:b83e:1f1a:a8ca References: User-Agent: G2/1.0 MIME-Version: 1.0 Message-ID: <7fe2291a-bc12-4708-85aa-0ffbdc25b2bfn@googlegroups.com> Subject: Re: Quick inverse square root From: Matt Borchers Injection-Date: Sun, 03 Jan 2021 22:31:16 +0000 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable Xref: reader02.eternal-september.org comp.lang.ada:61020 List-Id: 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:=20 > > I'm sure many of you have seen the Fast Inverse SquareRoot algorithm fr= om the open source Quake III engine. I just encountered it a few days ago. = Here it is, a bit reduced, from the original source:=20 > >=20 > > //C code from Quake III engine=20 > > float Q_rsqrt( float number )=20 > > {=20 > > long i;=20 > > float x2, y;=20 > > const float threehalfs =3D 1.5F;=20 > >=20 > > x2 =3D number * 0.5F;=20 > > y =3D number;=20 > > i =3D *(long *) &y;=20 > > i =3D 0x5f3759df - ( i >> 1 );=20 > > y =3D *(float *) &i;=20 > > y =3D y * (threehalfs - (x2 * y * y));=20 > > return y;=20 > > }=20 > >=20 > > It is interesting how much clearer the Ada code version is:=20 > >=20 > > with Interfaces; use Interfaces;=20 > > function QUICK_INVERSE_SQRT( a : FLOAT ) return FLOAT is=20 > > y : FLOAT :=3D a;=20 > > i : UNSIGNED_32;=20 > > for i'Address use y'Address;=20 > > begin=20 > > i :=3D 16#5F3759DF# - shift_right( i, 1 );=20 > > return y * (1.5 - (0.5 * a * y * y));=20 > > 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:=20 > >=20 > > 3/2 * 2**23 * (127 - mu) where mu is a constant close to 0.043.=20 > >=20 > > 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 calculate= d from the same formula but adjusted for the Long_Float bit layout.=20 > >=20 > > 3/2 * 2**52 * (1023 - mu) where mu is the identical constant as used fo= r Float case.=20 > >=20 > > 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=20 > calculated from the exponent must be numerically same for both float=20 > (IEEE 754 single-precision floating-point) and double (IEEE 754=20 > single-precision floating-point). Which is apparently not. You should=20 > get the exponent multiplied by the same power of 2 as for float. For=20 > double, I make a wild guess, you should replace right shift by 1 with=20 > right shift by 30 =3D 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 we= ll as the hex constant.=20 > >=20 > > generic=20 > > type F is digits <>;=20 > > type U is mod <>;=20 > > magic : U;=20 > > function G_Q_INV_SQRT( a : F ) return F;=20 > >=20 > > I write the body like this:=20 > >=20 > > function G_Q_INV_SQRT( a : F ) return F is=20 > > y : F :=3D a;=20 > > i : U;=20 > > for i'Address use y'Address;=20 > > begin=20 > > i :=3D magic - shift_right( i, 1 );=20 > > return y * (1.5 - (0.5 * a * y * y));=20 > > end G_Q_INV_SQRT;=20 > >=20 > > function QUICK_INVERSE_SQRT is=20 > > new G_Q_INV_SQRT( FLOAT, UNSIGNED_32, 16#5F3759DF# );=20 > >=20 > > This won't compile because the type U is not valid for the call to "shi= ft_right". How do I overcome this obstacle?=20 > >=20 > > 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 gen= eric at compile time? > As Jeffrey already suggested, use S'Exponent to extract the exponent.=20 >=20 > General notes.=20 >=20 > 1. C code relies on float being IEEE 754 single-precision floating-point= =20 > number with endianness opposite to integer endianness numbers. The=20 > exponent must land in the integer's MSB. This is clearly non-portable.=20 >=20 > 2. The approximation is very crude. I am too lazy to estimate its=20 > precision within the intended range, which is what? [0, 1]?=20 >=20 > 3. Ergo, making it generic has no sense.=20 >=20 > 4. If you port it to Ada. Add assertions validating endianness and=20 > floating-point format.=20 >=20 > --=20 > Regards,=20 > Dmitry A. Kazakov=20 > 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 :=3D 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 in= stantiate the generic and hope to find a way to eliminate it. I have not l= ooked 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 par= ameter. 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 proj= ect. 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 th= at 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 co= mpiler required me to instantiate the generic with different names and then= use renames for the function in the package specification.=20 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 :=3D U(1.5 * 2.0**(F'Machine_Mantissa - 1) * (F(= F'Machine_Emax - 1) - 0.043)); 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_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;