fortran 95 rounding up on it's own

Question

I decided to learn the fortran95 language (reason why is not important). However being a beginner I ran into a weird problem I really can't explain, therefor I need help.

I have the insertion sort algorithm:

subroutine insertion_sort_REAL4(array, array_len)
   implicit none
!parameners
   integer :: array_len
   real (kind=4), dimension(array_len) :: array 
!variables
   integer :: i,key,hole_pos
   do i = 0,array_len
      key = array(i)
      hole_pos = i;
      do while ((hole_pos > 0.0) .and. (key < array(hole_pos - 1)))
         array(hole_pos) = array(hole_pos - 1)
         hole_pos = hole_pos - 1
      end do
      array(hole_pos) = key
   end do
   return
end   

And there is the main program (excerpt):

real (kind = 4), dimension(3) :: x
x(1) = 3.1
x(2) = 4.3
x(3) = 5.4
write(*,*) 'Array = ',x
call insertion_sort_REAL4(x,3)
write(*,*) 'Array = ',x  

The first write statement prints out

Array =    3.09999990       4.30000019       5.40000010 

Why have the numbers been slightly changed? Does fortran95 not use the IEEE754 standard by default?

But let's say I can live with the slight alteration; the second write statements prints out

Array =    3.00000000       4.00000000       5.00000000  

Why have the numbers been rounded up? It's really bugging me and formatting the 'write' statement isn't doing any good and the Google searches haven't really helped. I guess there is not that many things on the internet about fortran as it is of C. I am a decent C programmer so any parallels to it are appreciated. Thanks you for your help!

Answer 1

A decimal number such as "3.1" will likely not have an exact representation in a binary number of finite length. The source code statement x(1) = 3.1 causes the computer to convert that decimal number into binary and store it. The statement write (*, *) x(1) causes the computer to fetch this binary value and convert it to decimal. Because "3.1" could not be exactly represented in finite-length binary, the conversion to decimal does not precisely recover "3.1". This explains the output of "3.09999990". This is not Fortran specific, but general to finite precision floating point arithmetic.

As to the other problem, key is declared integer in the sort subroutine and is thus rounding the reals to integers. When I compiled your program with full compiler warnings turned on, gfortran notified me of this.

If you you gfortran, try the following compiler options: -O2 -fimplicit-none -Wall -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wimplicit-interface -Wunused-parameter -fwhole-file -fcheck=all -std=f2008 -pedantic -fbacktrace. You will also find that your program has a subscript error.

Answer 2

For the first part: it does use IEEE754, and that is the cause why the numbers are "changed".

The What Every Computer Scientist Should Know About Floating-Point Arithmetic article is a must read to understand how this is working, and there are good IEEE754 calculators too...

So 3.1 was never exactly 3.1, but

3.0999999046325684

in the first place.

As for the second part: they are not rounded up but converted to integers, but I'm not into Fortran, so I guess something is declared as int in the insertion_sort_REAL4 routine, that causes the numbers to be converted to integers.