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 autolearn=ham autolearn_force=no version=3.4.5-pre1 Path: eternal-september.org!reader02.eternal-september.org!aioe.org!5WHqCw2XxjHb2npjM9GYbw.user.gioia.aioe.org.POSTED!not-for-mail From: "Dmitry A. Kazakov" Newsgroups: comp.lang.ada Subject: Re: Quick inverse square root Date: Sun, 3 Jan 2021 11:58:38 +0100 Organization: Aioe.org NNTP Server Message-ID: References: NNTP-Posting-Host: 5WHqCw2XxjHb2npjM9GYbw.user.gioia.aioe.org Mime-Version: 1.0 Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit X-Complaints-To: abuse@aioe.org User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:78.0) Gecko/20100101 Thunderbird/78.6.0 Content-Language: en-US X-Notice: Filtered by postfilter v. 0.9.2 Xref: reader02.eternal-september.org comp.lang.ada:61019 List-Id: 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