Re: Fortran compatibility

[Another Loco y ya <0xccaaffee@gmail.com>]

Dear Bill, and friends here, have you ever seen it? this ( C? ) code:

    float InvSqrt (float x){
     float xhalf = 0.5f*x;
     int i = *(int*)&x;
     i = 0x5f3759df - (i>>1);
     x = *(float*)&i;
     x = x*(1.5f - xhalf*x*x);
     return x;
    }

At some point of history, it is supposed to have been used widely for
calculating the 1 / sqrt(x).

i still wonder how it works, and inspired by the FORTRAN question you
just solved, i thought this would be worthy of some contemplation...
perhaps later. (At least me) Not hurried.

:-) Merci

2020-09-18 13:49 GMT-04:00, Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>:
> On Fri, Sep 18, 2020 at 04:39:47PM +0200, Bill Allombert wrote:
>> PROGRAM prog
>>   use ISO_C_BINDING, only : C_PTR, C_DOUBLE
>>   use PARI
>>   implicit none
>>   real(kind=C_DOUBLE) :: r      = 1e36
>>   type(C_PTR)         :: p
>>   CALL pari_init(10000000_8,2_8)
>>   p = gmul(gprec(dbltor(r),1000_8),glog(stoi(10000_8),20_8))
>>   !p = gmod(p, Pi2n(1_8,20_8))
>>   r = rtodbl(p)
>>   CALL pari_close()
>>   PRINT '(a,f0.9)','1e36*log(10000)%(2*Pi)  = ', r
>> END PROGRAM prog
>
> Sorry my mailer sent the wrong file, the correct version is
>
> PROGRAM prog
>   use ISO_C_BINDING, only : C_PTR, C_DOUBLE
>   use PARI
>   implicit none
>   real(kind=C_DOUBLE) :: r      = 1e36
>   type(C_PTR)         :: p
>   integer(kind=C_LONG) :: prec = 20 ! 18 words
>   CALL pari_init(10000000_8,2_8)
>   p = glog(stoi(10000_8),prec)
>   p = gmod(p, Pi2n(1_8, prec))
>   r = rtodbl(p)
>   CALL pari_close()
>   PRINT '(a,f0.9)','log(10000)%(2*Pi)  = ', r
> END PROGRAM prog
>
> Cheers,
> Bill
>
>