fortran - Fortran 中非常小的数字变成负数

Question

我正在用 Fortran90 编写一个程序，它计算从 i=1 到 i=n 的总和，其中给出了 n。总和为 sum_{i=1}^{i=n}1/(i*(i+1)*(i+2))。这个总和收敛到 0.25。这是代码:

PROGRAM main

INTEGER n(4)
DOUBLE PRECISION s(4)
INTEGER i

  
OPEN(11,FILE='input')
OPEN(12,FILE='output')
DO i=1,4
  READ(11,*) n(i)
END DO

PRINT*,n
CALL suma(n,s)
PRINT*, s

END

 
SUBROUTINE suma(n,s)
INTEGER n(4),j,k
DOUBLE PRECISION s(4),add
s=0 
DO k=1,4
  DO j=1,n(k)
    add=1./(j*(j+1)*(j+2))
    s(k)=s(k)+add 
  END DO 
END DO
  
END SUBROUTINE

输入

178   
1586     
18232 
142705

output 文件现在是空的，我需要对它进行编码。我只是打印结果，它们是: 0.249984481688 0.249999400246 0.248687836759 0.247565846142

问题出在变量 add 上。当j较大时，add变为负数，求和收敛不好。我该如何解决？

Answer 1

问题是整数溢出。 142705142706142707 是一个对于 4 字节整数来说太大的数字。

然后发生的是数字溢出并循环回负数。

正如@albert 在他的评论中所说，一种解决方案是在每一步都将其转换为 double :((1.d0/j)/(j+1) )/(j+2)。这样，它就是用浮点值进行计算。

另一种选择是使用 8 字节整数:

integer, parameter :: int64 = selected_int_kind(17)
integer(kind=int64) :: j

不过，您在计算时应该非常小心。更好并不总是更好。我建议您查看计算机如何执行浮点运算，以及这会产生什么问题。参见示例 here on wikipedia .