[Bug fortran/54159] Fortran quad precision rounding seemingly nonsensical

[dominiq at lps dot ens.fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=54159

Dominique d'Humieres <dominiq at lps dot ens.fr> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |RESOLVED
         Resolution|                            |INVALID

--- Comment #3 from Dominique d'Humieres <dominiq at lps dot ens.fr> 2012-08-02 11:35:32 UTC ---
What you see is the "round to even" effect as shown with the simpler following
test:

program quad_test
   implicit none

   integer, parameter :: dp = kind(0d0)
   integer, parameter :: qp = selected_real_kind(32)

   real(qp) :: y

   y = 1.0_qp+0.5*real(epsilon(1.0_dp), kind=qp)
   print *, y, nearest(y,1.0_qp)
   print *, real(y, kind=dp), real(nearest(y,1.0_qp), kind=dp)
end program quad_test

which outputs

   1.00000000000000011102230246251565404        
1.00000000000000011102230246251565423      
   1.0000000000000000        1.0000000000000002     

i.e., y is rounded to 1.0_dp (the closest even binary representation) while
nearest(y,1.0_qp) is rounded to the closest double: nearest(1.0_dp,1.0_dp). For
your numbers the hexa representations are:

400E4A0CCB6E8DB58800000000000001
400E4A0CCB6E8DB58800000000000000

40E4A0CCB6E8DB59
40E4A0CCB6E8DB58

and if I am not mistaken, the (default) rounding to "even" are correct.

Your comments "Two quad prec numbers that differ by last 3 digits" and
"truncated output" make me suspect that you are not fully aware of the detailed
behavior of floating-point representations and their default rounding
properties.