I do hydrodynamics simulations with Fortran and recently I met with this issue:
I have a single-precision array b
of length n=512**3
, and I wish to compute the root-mean-squre of its elements. All I do is
s=0
do i=1,n
s=s+b(i)**2
enddo
and then compute sqrt(s/n)
.
I found that it makes about $30\%$ difference in the result depending on whether I have declared s
as single-precision (real
) or double-precision (real*8
). My question is which should be the correct way to do?
My naive understanding is that, when I use real*8 s
, the code will first convert b(i)
into double precision and then do the computation. In that case, single-precision values in b
will be padded with zeros in base 2, which makes them slightly change in base 10. Even so, I thought it was still accurate in base 2 and there should be no harm of promoting them to double precision before the computation. On the other hand, if I stick to single precision then each time I take the square b(i)**2
, there is some tiny error introduced. Based on these arguments I though real*8
should be the right choice.
If my thoughts above were correct, does it mean one need to go for quadruple precision when dealing with double-precision data?
I'm pretty new to numerical computations and Fortran, so any answers/comments will be helpful!
binary32
tobinary64
is exact. An error of this kind may have been introduced earlier, in conversions from strings or literal constants typed inside your source code. $\endgroup$