print *, 1.0+epsilon(1.0) ! should identify significant digits
end
No need to repost an old incorrect answer. Your code just prints the number 1.00000012, which can be represented as 1.0000001 without loss. So it only needs 8 digits, not 9. But there are numbers that really need 9 digits, like the one I gave in my answer.
@certik suggested print *, x, y for checking the values of x and y. Different compilers choose different numbers of digits to print (the Fortran standard allows that), and so I recommend using an explicit format such as '(F0.8)' or even '(F0.9)' etc. when one is investigating the last digit.
This is sometimes a subtle question. In fortran the precision intrinsic returns an integer for each real kind. Precision(1.0) return 6, which means that every six-decimal-digit value (within the exponent range) can be represented correctly with rounding. It does not mean that those decimal numbers are represented exactly, it means that they can be distinguished correctly with rounding from their adjacent decimal numbers. In contrast, epsilon(1.0) returns the difference between 1.0 and the next largest binary floating point number. For ieee real32, that value is about 1.2e-7. The precision is 6 because the epsilon is just a little too large for it to be 7. The other feature that has been mentioned in this discussion is the number of decimal digits that are necessary to distinguish, with rounding, between adjacent binary numbers. That number is typically about 3 larger than the precision.
I am not sure if I understand it correctly. But I think it is compiler dependent too.
For this code on windows 10:
program sp
implicit none
real(4) :: v = 1.123456789
print*, v
end program sp
with gfortran version (13.2.0)
gfortran main.f90 -o main
i get
1.12345684
with intel
Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2024.0.0 Build 20231017
Copyright (C) 1985-2023 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 14.32.31332.0
Copyright (C) Microsoft Corporation. All rights reserved.
I get
1.123457
This is just due to the list directed i/o. A compiler is free to print as many digits as it wants, with arbitrary spacing between values, and with arbitrary line feeds. So even if the binary representations are the same (which you could verify with a hex output format), the list directed outputs can look very different.
2 Likes
That is why I said what I did yesterday but @Shahid seems not to have noticed it. This program reveals the danger of list-directed (format *) i/o:
implicit none
integer, parameter:: dp = kind(1d0)
real(dp) :: bigin = huge(1.0_dp), bigout
character:: charbig*32, msg*200
integer :: io
write(charbig,*) bigin
print "(A,ES23.15E3,2A)",'bigin =',bigin,' charbig = ',charbig
print "(A,ES24.16E3,2A)",'bigin =',bigin,' charbig = ',charbig
read (charbig,*,iostat=io,iomsg=msg) bigout
if(io==0) then
print *,'bigout=',bigout
else
print *,trim(msg)
end if
end program
Compiling it with ifx gave this at run time:
bigin = 1.797693134862316E+308 charbig = 1.797693134862316E+308
bigin = 1.7976931348623157E+308 charbig = 1.797693134862316E+308
input conversion error, unit -5, file Internal List-Directed Read
That is because the ifx internal write, which is standard-conforming, produced a value of charbig that was bigger than huge(1.0_dp). Gfortran used one more significant digit and ran happily. That is also standard-conforming.
1 Like
I think this is a somewhat different issue than simply rounding up or down correctly. If the value of huge(1.0_wp) is rounded up at any decimal place, then that value cannot be represented in the floating point representation of that kind. I don’t think it matters if it is list-directed i/o or formatted i/o, the value will be too large to convert. Consider the following program:
program round
use, intrinsic :: iso_fortran_env, only: wp=>real64, int64
implicit none
call readx('1.79769313486231570814527423731704+308')
call readx('1.79769313486231570814527423732+308')
call readx('1.79769313486231570814527424+308')
call readx('1.7976931348623157081453+308')
call readx('1.79769313486231570815+308')
call readx('1.79769313486231571+308')
call readx('1.797693134862316E+308')
call readx('1.7976931349E+308')
call readx('1.797693135E+308')
call readx('1.7977E+308')
call readx('1.798E+308')
call readx('1.8E+308')
call readx('2.E+308')
contains
subroutine readx( cx )
character(*), intent(in) :: cx
real(wp) :: x
integer :: istat
read(cx,*,iostat=istat) x
if ( istat .ne. 0 ) then
write(*,'("istat=",i0," cx=",a)') istat, cx
else
write(*,'("x=",z16.16," cx=",a)') transfer(x,1_int64), cx
endif
return
end subroutine readx
end program round
$ nagfor round.f90 && a.out
NAG Fortran Compiler Release 7.2(Shin-Urayasu) Build 7200
[NAG Fortran Compiler normal termination]
x=7FEFFFFFFFFFFFFF cx=1.79769313486231570814527423731704+308
x=7FEFFFFFFFFFFFFF cx=1.79769313486231570814527423732+308
x=7FEFFFFFFFFFFFFF cx=1.79769313486231570814527424+308
x=7FEFFFFFFFFFFFFF cx=1.7976931348623157081453+308
x=7FEFFFFFFFFFFFFF cx=1.79769313486231570815+308
x=7FEFFFFFFFFFFFFF cx=1.79769313486231571+308
istat=215 cx=1.797693134862316E+308
istat=215 cx=1.7976931349E+308
istat=215 cx=1.797693135E+308
istat=215 cx=1.7977E+308
istat=215 cx=1.798E+308
istat=215 cx=1.8E+308
istat=215 cx=2.E+308
That first value is huge(1.0_wp) with a lot of decimals. I just did the rounding for the other values from that by hand, but I think they are all correct. All of those rounded up values will exceed huge, so they could all generate a read i/o error. The nagfor output is shown above, and it does produce i/o errors for 16 decimal digits and less, but it does not error for the cases with more digits. That behavior cannot be explained without knowing exactly how the character strings are converted to the floating point values, but one can imagine that the input digits might ignored beyond the point at which they could affect the last bit in the binary value. Or maybe the conversion is done from right to left (least significant to most significant digits), which means that the least significant digits are numerically ignored during the conversion process. Or maybe the conversion is done in real128 arithmetic, and the errors are detected when converted to real64. There are several possible approaches that would result in this output.