Output [with WRITE(*,*)] of floating point numbers does not comply with the IEEE-754 Standard

[Hans_Peschke]

Hi,

according to the IEEE Standard for Floating-Point Arithmetic [1] (especially
clause 5.12 paragraph 2) a conversion from a binary format to external decimal
format and back (assumed that the fp number is representable in all involved
formats and the rounding mode is nearest) has to recover the original value
(sign, mantissa and exponent).

Imho this should be especially the case for the 'default output statement.

But this is not the case, if the WRITE statement is used with '*' as the
format specification and compiled with ifort (Linux 32-Bit). There is missing
at least one extra digit to comply with the standard. This is true for single,
double and quadrupel precision.

for example:

[fortran]real(kind(0.0D0)) :: x, y

x = tiny(1.0D0) ! = ulp(1.0)

write(*,*) x ! -> 2.225073858507201E-308
! correct value:  2.2250738585072014E-308 ! with all significand decimal digits

! read in the same character string and convert to double precision fp number y
read(*,*) y  ! <- 2.225073858507201E-308

! this results in
! mant(x) == 0.50000000000000000E-000
! mant(y) == 0.22250738585072009E-307[/fortran]

The following program enables you to easily test the behaviour with the available
fp precisions (single, double and quadrupel for ifort and single, double and
extended for g95 and gfortran).

[fortran]program ioconv
implicit none

! kinds for real floating point precision numbers
integer, parameter :: sp = kind(0.0)   &                ! single
                    , dp = kind(0.0D0) &                ! double
                    , ep = selected_real_kind(18, 1000) ! extended (g95 & gfortran)
integer, parameter :: qp = kind(0.0Q0)                  ! quadrupel (ifort only)

! format strings for exact output (default)    s0.dd...ddEse...e
character(*), parameter :: sfmt = '(E16.9  E2)'   !  9 essential decimal digits
character(*), parameter :: dfmt = '(E25.17 E3)'   ! 17 essential decimal digits
character(*), parameter :: efmt = '(E30.21 E4)'   ! 21 essential decimal digits
character(*), parameter :: qfmt = '(E45.36 E4)'   ! 36 essential decimal digits

! format strings for exact output (scientific) sd.dd...ddEse...e
character(*), parameter :: sfmts = '(ES15.8  E2)' !  9 essential decimal digits
character(*), parameter :: dfmts = '(ES24.16 E3)' ! 17 essential decimal digits
character(*), parameter :: efmts = '(ES29.20 E4)' ! 21 essential decimal digits
character(*), parameter :: qfmts = '(ES44.35 E4)' ! 36 essential decimal digits

! please CHOOSE appropiate and matching values from constants above
integer, parameter      :: prec     = dp
character(*), parameter :: precfmts = dfmts
character(*), parameter :: precfmt  = dfmt

! some digits of pi
real(prec), parameter :: pi = 3.1415926535897932384626433832795028841971693993751_prec

real(prec) :: z

!read(*,*) z
z = tiny(0.0_prec)
!z=epsilon(0.0_prec)
!z = pi

call iotest(z)

contains

  subroutine iotest(x)
    real(prec) :: x, y
    
    write(*,*) 'x = '
                    WRITE(*,*)            x
    call onechar(); WRITE(*,fmt=precfmts) x
    call onechar(); WRITE(*,fmt=precfmt)  x

    write(*,fmt='(A)') '   ^    ^    ^    ^    ^    ^    ^    ^    ^'
    write(*,fmt='(A)') '   0    5   10   15   20   25   30   35   40'

    do
      write(*,*)
      write(*,*) 'Please copy and paste one of the above fp numbers to check, ', &
                 'if the printed digits are enough to represent the correct ',   &
                 'value of x. (stop with 1.0)'

      write(*,'("y := ")', advance='no'); read(*,*) y

      write(*,*) '(x == y) = ', x == y
      write(*,*) 'mantissas and error: (fraction(x), fraction(y), x-y):'
      write(*,fmt=precfmt) fraction(x)
      write(*,fmt=precfmt) fraction(y)
      write(*,fmt=precfmts) x - y 
      if (y == 1.0_prec) return
    end do
  end subroutine iotest

  ! print only one char, only used for output allignment
  subroutine onechar
    write(*,fmt='(" ")', advance='no')
  end subroutine onechar
end program ioconv[/fortran]

Best regards

Hans

Links:
[1] http://grouper.ieee.org/groups/754/

[TimP]

According to my reading of IEEE754, the assurance of full accuracy in this situation applies only when you specify the necessary number of digits (17 for double precision), and in general only for magnitudes in a limited range such as < 0.1d0 < |x| < 1/epsilon(x). IEEE754 only requires there be a way to produce sufficient digits, not that it be a default. There is no assurance with list-directed format; this is generally considered to be the more popular choice, as th 17th digit is only seldom significant.