为什么同一程序的c和fortran版本会产生不同的结果？

编程问答 2022-06-03

问题描述

我使用gcc（Ubuntu 9.3.0-17ubuntu1〜20.04）9.3.0

c代码是：

// 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代码：

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

program little

implicit none

integer,parameter  :: n = 800
real                :: a(n,n),b(n,c(n,n)
real                :: t
integer             :: i,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
        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,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

c程序打印：1362136192.000000

和Fortran程序打印：1362137216.0000

如果我不将每个元素本身相乘，正如我在代码的注释中指出的那样，则对于程序的两个版本，我将获得相同的值：

c图元：639200.000000

Fortran程序：639200.0000

为什么当我使用乘法时，c和Fortran代码会产生不同的结果？是否必须使用不同的实数实现？

解决方法

差异是由于求值顺序与浮点型精度的限制所致。

如果将Fortran版本更改为

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

即在a和b的两边加上括号，两种语言的结果相同。由于使用了+=赋值运算符，C版本已经使用了此评估顺序。

如评论中所述，这是在可用精度限制下的预期行为。

,

当我编写程序的Ada版本时，我发现必须将小数精度降低到6个小数才能获得Fortran答案。

Ada版本为：

和Ada.Text_IO;使用Ada.Text_Io; 

procedure Main is
   type Real is digits 6;
   package Real_IO is new Ada.Text_IO.Float_IO(Real);
   use Real_IO;
   
   subtype Index is integer range 1..800;
   type Matrix is array(Index,Index) of Real;
   
   A : Matrix;
   B : Matrix;
   C : Matrix;
   T : Real := 0.0;
begin
   for I in Index loop
      for J in Index loop
         A(I,J) := Real(I - 1 + J - 1);
         B(I,J) := Real(I - 1 - (J - 1));
      end loop;
   end loop;
   
   for I in Index loop
      for J in Index loop
         T := 0.0;
         for K in Index loop
            T := T + A(I,K) * A(I,K) + B(K,J) *B(K,J);
         end loop;
         C(I,J) := T;
      end loop;
   end loop;
   
   Put(Item => C(Index'Last,Index'Last),Exp => 0,Aft => 4);
   New_Line;
end Main;

定义类型Real的行定义浮点类型的精度：

type Real is digits 6;

使用六位精度产生的值是

1362137216.0000

使用更高精度的浮点类型会得出该值

1362135200.0000
cfortrangfortranprecision