提问人:Antonio Serrano 提问时间:11/12/2020 更新时间:11/12/2020 访问量:241
为什么同一程序的 c 和 fortran 版本会产生不同的结果?
Why the c and fortran versions of this same program produce different results?
问:
我用的是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,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
c 程序打印:1362136192.000000
Fortran 程序打印:1362137216.0000
如果我不将每个元素本身相乘,正如我在代码的注释中所述,我会得到两个版本的程序相同的值:
c prigram: 639200.000000
Fortran 计划:639200.0000
为什么当我使用乘法时,c 和 Fortran 代码会产生不同的结果?它是否必须与实数的不同实现有关?
答:
6赞
rtoijala
11/12/2020
#1
这种差异是由于计算顺序与浮点类型的有限精度相结合所致。
如果将 Fortran 版本更改为
t = t + (a(i,k) * a(i,k) + b(k,j) * b(k,j))
即在术语 和 周围添加括号,您将获得两种语言的相同结果。由于使用了赋值运算符,C 版本已经使用了此计算顺序。a
b
+=
如注释中所述,这是在可用精度限制下的预期行为。
评论
1赞
Eric Postpischil
11/12/2020
或者,在我的 C 实现中,更改 C 版本中的行以获得与 FORTRAN 版本相同的结果(Apple Clang 11,带有 )。我怀疑可能有多种评估方法可以得到相同的结果。t += b[k][j] * b[k][j] + a[i][k] * a[i][k];
-O3
1赞
John Alexiou
11/12/2020
请注意,两个不同量级的加法和最重要的减法会引入明显的舍入误差。
0赞
Jim Rogers
11/12/2020
#2
当我编写该程序的 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
评论
0赞
Vladimir F Героям слава
11/12/2020
那么,这个问题的答案是什么?
0赞
Jim Rogers
11/12/2020
答案是,差异可能由多种因素引起,包括由于不同的计算顺序而导致的不同舍入误差,以及由于精度降低而导致的舍入误差,浮点数。
1赞
Vladimir F Героям слава
11/12/2020
但是问题中的 Fortran 和 C 都使用单精度。证明它在不同语言中的双精度会有所不同有什么意义?很明显,在任何语言的双精度下,它都会有所不同。
0赞
Jim Rogers
11/12/2020
我使用了 GNAT Ada 编译器,它是 GNU 工具链的一部分。它应该产生与GNU C++编译器相同的结果。两种语言共享相同的编译器后端。
评论
float
double
float