在 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*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。

对我来说，在重载运算符周围胡思乱想是个坏主意。如果确实有太多代码需要手动修复，请查看您的编译器是否可以选择将默认文字类型更改为双精度。