How come when I assign a simple value (or even an expression) to a value I sometimes get less digits of precision than I expect?
In the following example program the DOUBLEPRECISION value VALUE2 only contains the same number of correct digits from the constant that was assigned to it as the REAL value VALUE1. VALUE3 contains significantly more...
module M_prec
integer,parameter :: sp = selected_real_kind(6)
integer,parameter :: dp = selected_real_kind(15)
end module M_prec
program main
use M_prec
real :: value1=123.45678901234567890123
doubleprecision :: value2=123.45678901234567890123
doubleprecision :: value3=123.45678901234567890123d0
real(kind=dp) :: value4=123.45678901234567890123_dp
write(*,*)value1
write(*,*)value2
write(*,*)value3
write(*,*)value4
end
123.456787
123.45678710937500
^^^^^^^^ <== Not the values you might expect
123.45678901234568
123.45678901234568
A general principle of Fortran is that the type of the RHS (Right Hand Side) of an assignment does not depend on the LHS (Left Hand Side). Once you understand this rule, a lot of things fall into place.
You must make sure the constant on the RHS is of a KIND that can hold the number of digits desired, regardless of what type the LHS is. "123.45678901234567890123" is a REAL expression and is evaluated first. Then it is assigned to the value on the LHS, which in this case promotes the REAL value to DOUBLEPRECISION. Adding the "d0" suffix or specifying and using a KIND sufficiently big enough to give you the accuracy you desire is required.
The best explanation I have seen for this is by Dick Hendrickson on the newsgroup "comp.lang.fortran" :
I think the REAL reason (ho, ho, ho) is that your example is too simple. In Fortran almost everything can be an expression, rather than a single term. (heck, even 29.5... is an [degenerate] expression in the syntax rules) And there is no good way to push the type of the left hand side into the right hand side expressions that won't surprise someone. In something like
double_prec_var = N/4 + MIN(M/3, 10) + 7 * user_func (6, 3.14)everything goes wrong if you pull the constants up to double precision before evaluating the RHS. Especially if userfunc is generic or if it returns a type that causes the asterisk in "7 * userfunc (6, 3.14)" to be a user defined function (probably also generic).
Rather than have two sets of rules (either for "expressions" that are constants or for expressions that appear in declaration statements) Fortran chose one rule for all expressions and it confuses people on the simple case. It's a pity and, IMO, compilers should be more aggressive about mismatched precision with simple constant expressions.
Note that compilers are free to produce a warning when all the digits of a constant are not stored. As an example, if you use the -Wconversion-extra switch on gfortran:
+ gfortran -Wconversion-extra xxx.f90
xxx.f90:8:52:
real :: value1=123.45678901234567890123
1
Warning: Non-significant digits in 'REAL(4)' number at (1), maybe incorrect KIND [-Wconversion-extra]
xxx.f90:9:28:
doubleprecision :: value2=123.45678901234567890123
1
Warning: Conversion from 'REAL(4)' to 'REAL(8)' at (1) [-Wconversion-extra]
xxx.f90:9:52:
doubleprecision :: value2=123.45678901234567890123
1
Warning: Non-significant digits in 'REAL(4)' number at (1), maybe incorrect KIND [-Wconversion-extra]
xxx.f90:10:54:
doubleprecision :: value3=123.45678901234567890123d0
1
Warning: Non-significant digits in 'REAL(8)' number at (1), maybe incorrect KIND [-Wconversion-extra]
xxx.f90:11:55:
real(kind=dp) :: value4=123.45678901234567890123_dp
1
Warning: Non-significant digits in 'REAL(8)' number at (1), maybe incorrect KIND [-Wconversion-extra]