Fortran 中非常少的数字变为负数答案

【问题标题】：Very small number turns negative in FortranFortran 中非常少的数字变为负数
【发布时间】：2020-11-12 10:22:08
【问题描述】：

我正在 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

输入

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

问题在于变量add。当j 较大时，add 变为负数，并且总和不能很好地收敛。我该如何解决？

【问题讨论】：

@KenY-N 不是在READ(11,*) n(i) 中初始化的n？
文件input中的值是什么？请将文件添加到问题中。结果是什么（PRINT*,n 和PRINT*, s 的结果）？也请添加到问题中。
是j*(j+1)*(j+2)溢出整数变成负数的问题吗？如果n(k)大约大于最大整数值的立方根，那么可能存在溢出，导致负数？
根本不是一个好的解决方案，但可以写：add=((1./j)/(j+1))/(j+2)（免责声明，未测试）。
由于 KenY-N 给出的原因而失败，1./j*(j+1)*(j+2) 术语中存在整数溢出，因此环绕为负数。请注意，您的计算是在real 中执行的，而不是在double precision 中执行的。更好的是add=((1.d0/j)/(j+1))/(j+2)

标签： fortran

【解决方案1】：

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

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

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

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

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

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

【讨论】：

我没有看到 int64 位使用的真正优势，(j*(j+1)*(j+2) 中的这种溢出不会也很快（可能在 2^21 左右，只有 200 万）？
同意@albert 除非它对性能至关重要，否则我想我会先除以 j，然后除以 (j + 1)，然后除以 (j + 2)
如果您想避免额外的除法，Fortran 提供了转换内在函数。试试1/(real(j) * real(j+1) * real(j+2))。更好的选择可能是递归计算加数：a(1) = 1 / real(1 * 2 * 3) with a(j+1) = j * a(j) / (j + 3)`
如果你真的想避免除法，我会预先计算逆因子直到 n(k) 的最大值，并在内部循环中使用乘法

【解决方案2】：

这可能是实现您想要的更好的方法。我确实删除了 IO。程序的输出是

% gfortran -o z a.f90 && ./z
   178 0.249984481688392
  1586 0.249999801599584
 18232 0.249999998496064
142705 0.249999999975453

program main

  implicit none  ! Never write a program without this statement

  integer, parameter :: knd = kind(1.d0) ! double precision kind

  integer n(4)
  real(knd) s(4)
  integer i

  n = [178, 1586, 18232, 142705]
  call suma(n, s)
  do i = 1, 4
    print '(I6,F18.15)', n(i), s(i)
  end do

  contains
     !
     ! Recursively, sum a(j+1) = j * a(j) / (j + 1)
     !
     subroutine suma(n, s)
        integer, intent(in) :: n(4)
        real(knd), intent(out) :: s(4)
        real(knd) aj
        integer j, k
        s = 0 
        do k = 1, 4
           aj = 1 / real(1 * 2 * 3, knd)    ! a(1)
           do j = 1, n(k)
              s(k) = s(k) + aj
              aj = j * aj / (j + 3)
           end do 
        end do
     end subroutine

end program main

【讨论】：