Why the c and fortran versions of this same program produce different results?

Question

I use gcc (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0

The c code is:

// Compile with:
//  gcc -o little_c little.c
#include <stdio.h>  // printf

void main(void) {
    int n = 800;
    float a[n][n], b[n][n], c[n][n];
    int i, j, k;

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            a[i][j] = (float) (i+j);
            b[i][j] = (float) (i-j);
        }
    }

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            float t = (float) 0.0;
            for (k = 0; k < n; k++)
                t += a[i][k] * a[i][k] + b[k][j] * b[k][j];
                //t += a[i][k] + b[k][j]; // If I comment the above line and uncomment this, the c and fortran reults are the same
            c[i][j] = t;
        }
    }
    printf("%f", c[n-1][n-1]); // prints the very last element

}

Fortran code:

! Compile with:
!  gfortran -o little_fort little.f90 

program little

implicit none

integer, parameter  :: n = 800
real                :: a(n,n), b(n,n), c(n,n)
real                :: t
integer             :: i, j, k  ! Counters

do i = 1, n
    do j = 1, n
        a(i,j) = real(i-1+j-1)      ! Minus one, for it to be like the c version
        b(i,j) = real(i-1-(j-1))    ! because in c, the index goes from 0 to n-1
    end do
end do

do i = 1, n
    do j = 1, n
        t = 0.0
        do k = 1, n
            t = t + a(i,k) * a(i,k) + b(k,j) * b(k,j)
            !t = t + a(i,k) + b(k,j) ! If I comment the above line and uncomment this, the c and fortran reults are the same
        end do
        c(i,j) = t
    end do
end do
    
write (*,"(F20.4)") c(n,n)  ! This is the same as c[n-1][n-1] in c


end program little

The c program prints: 1362136192.000000

and the Fortran program prints: 1362137216.0000

If I do not multiply each element by itself, as I state in the comments in the code, I get the same value for both versions of the program:

c prigram: 639200.000000

Fortran program: 639200.0000

Why when I use a multiplication the c and Fortran code produce different results?. Does it have to be with different implementations of the real numbers?

Answer 1

The difference is due to the order of evaluation combined with the limited precision of the floating point type.

If you change the Fortran version to

t = t + (a(i,k) * a(i,k) + b(k,j) * b(k,j))

i.e. add parenthesis around the terms with a and b, you get the same result for both languages. The C version already uses this order of evaluation due to the use of the += assignment operator.

As mentioned in the comments, this is expected behavior at the limits of the available precision.