Fortran中双精度除以整数

Question

现在我正在处理在我之前从事该项目的人的 Fortran 代码。代码比较大，这里就不提供了。在这个代码的很多地方都有双精度值除以整数。IE：

double precision :: a,c
integer :: b
a = 1.0d0
b = 3
c = a/b

这有时会导致比常规浮点数错误更严重的错误。我会为您复制粘贴其中一个这样的错误：3.00000032899483。对于他的目的来说，这是可以的，但对我来说，这个错误是可怕的。似乎fortran正在将实数除以整数，然后将数字放大到双精度，添加一些随机的东西。如果我会写，当然可以处理：

c = a/(real(b,8))

所以，我可以做到，但是搜索他的所有代码需要很长时间。这就是为什么，我尝试以下一种方式重载 (/) 运算符：

MODULE divisionoverload
!-----------------------------------------------------
   INTERFACE OPERATOR (/)
      MODULE PROCEDURE real_divoverinteger
   END INTERFACE OPERATOR ( / )
CONTAINS
!-----------------------------------------------------

DOUBLE PRECISION FUNCTION real_divoverinteger(a,b)
   IMPLICIT NONE
   DOUBLE PRECISION, INTENT (IN) :: a 
   INTEGER, INTENT(IN) :: b 

   real_divoverinteger = a/real(b,8)
END FUNCTION real_divoverinteger

END MODULE divisionoverload

但是 gfortran4.8 给我的错误清楚地表明这与内在除法运算符冲突：

 MODULE PROCEDURE real_divoverinteger
                                          1
Error: Operator interface at (1) conflicts with intrinsic interface

那么在不手动处理所有整数除法的情况下，我该怎么做才能解决这个问题呢？

使用打印格式生成此类数字的代码：

 hx = (gridx(Nx1+1)-gridx(Nx1))/modx)
 do i = 2,lenx
    sub%subgridx(i)=sub%subgridx(i-1) + hx
 end do
 !*WHERE*
 !hx,gridx[1d array], sub%subgridx[1d array] - double precision
 !modx, lenx - integer

打印代码

subroutine PrintGrid(grid,N)
    integer, intent (in) :: N       
    double precision, dimension(N), intent (in) :: grid
    integer :: i
    print"(15(f20.14))", (grid(i), i=1,N)
        print*, ""
end subroutine PrintGrid

---

解决方案： 我为这个问题道歉。我应该更仔细地阅读代码。我发现了一个错误。您遇到的问题的重现方式如下：

program main
   double precision :: a
   a = 1.0/3
   print*, a
end program main

所以基本上，他正在执行实数/整数除法并将这个值分配给 real*8。我用这样的代码改变了一切：

 program main
   double precision :: a
   a = 1.0d0/3
   print*, a
end program main

它奏效了！非常感谢您的帮助，如果没有您的回答，我将更多地处理这个问题。

Answer 1

我相信你的问题，这是最初的断言（如果我明白你也行）。将一个real由integer导致较低精度比一分real，是错误的。在计算操作之前，Fortran 将操作的两个操作数中最弱的类型（如果需要）提升为最强的类型。因此，自己做推广是没有必要的。

例如，请参阅以下代码：

program promotion
    implicit none
    integer, parameter :: dp = kind( 1.d0 )
    integer :: i
    real( kind = dp ) :: a

    i = 3
    a = 1
    print *, "division by a real:", a/real(i,dp), "division by an integer:", a/i
    print *, "are they equal?", a/real(i,dp) == a/i
end program promotion

结果的实际打印可能因一个编译器而异，但最终测试将始终评估为 .true.

~/tmp$ gfortran promotion.f90
~/tmp$ ./a.out 
 division by a real:  0.33333333333333331      division by an integer:  0.33333333333333331     
 are they equal? T