Re: Possible loss of accuracy

[Stephen Montgomery-Smith]

On 05/15/2013 11:59 AM, John W. Eaton wrote:
> On 05/15/2013 05:27 AM, Marco Caliari wrote:
>> Dear all,
>>
>> if I compute
>>
>> t=50
>> (1-2^(-t))^(2^t)
>>
>> I get (different versions, machines, ...)
>>
>> ans = 3.67879441128629e-01
>>
>> which is correct to the 10th digit (compared with www.wolframalpha.com).
>> I can stay with it, but Matlab is correct to the 15th digit. While
>> trying to understand, I wrote the following Fortran code
>>
>> program test
>> integer t
>> t = 50
>> write(6,'(1x,e24.16)') (1.d0-2.d0**(-t))**(2.d0**t)
>> write(6,'(1x,e24.16)') (1.d0-2.d0**(-50))**(2.d0**50)
>> stop
>> end
>>
>> If I compile with gfortran test.f (gfortran 4.3.3) and run, I get
>>
>> 0.3678794411286288E+00 (~ Octave value)
>> 0.3678794411714422E+00 (~ Matlab value)
>>
>> and if I compile with gfortran -O2 test.f and run, I get
>>
>> 0.3678794411714422E+00
>> 0.3678794411714422E+00
> 
> I don't know what the problem is, but with Octave 3.6.2 on my system,
> I get 0.367879441171442 for your expression.

Running octave under ubuntu linux I get 3.6787944112862880353e-01, but
with FreeBSD I get 3.6787944117144216749e-01.

1-2^(-50) should definitely be exactly representable in IEEE double
precision, and IEEE double precision should definitely be sufficient to
compute this to the required precision.  I think it is worth looking
inside glibc to see how exp and log are written.  I know that the code
used in FreeBSD is different, and probably is that written by Sun
Microsystems, and extreme care is taken to extract every last drop of
precision that IEEE can muster.