在 Fortran 中用整数除以双精度

Question

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

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/integer 除法并将这个值赋给 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

Answer 2

如前所述，您的标题声明不正确，整数除法会自动向上转换并且没问题。

你的第一个例子确实有问题：

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

这里我们有一个单精度1.。 3 被提升为单精度，因此您将 1/3 的单精度近似值分配给双精度 a。解决方法是使 1.0 1.d0。

我觉得绕过重载运算符是个坏主意。如果确实有太多代码需要手动修复，请查看您的编译器是否有将默认文字类型更改为双精度的选项。

在 Fortran 中用整数除以双精度

Division of double-precision by an integer in Fortran

precision

fortran

division

fortran90