使用不受舍入误差影响的 CEILING

Posted 2023-02-16

【中文标题】使用不受舍入误差影响的 CEILING

【英文标题】：Use CEILING Without the Effect of Rounding Error

【发布时间】：2016-10-11 03:09:47

【问题描述】：

我正在尝试使用内部函数“CEILING”，但舍入误差有时会导致难以获得我想要的结果。示例代码非常简单：

PROGRAM MAIN

  IMPLICIT NONE

  INTEGER, PARAMETER        ::      ppm_kind_double = KIND(1.0D0)
  REAL(ppm_kind_double)     ::      before,after,dx

  before = -0.112
  dx = 0.008
  after = CEILING(before/dx)

  WRITE(*,*) before, dx, before/dx, after

  END

我得到了结果：

我在代码中给'before'和'dx'的值只是为了演示。例如，对于之前/dx = -13.5 的那些，我想使用 CEILING 来获得 -13。但是对于我展示的图片，我实际上想要得到-14。我考虑过使用一些参数，例如

 IF(ABS(NINT(before/dx) - before/dx) < 0.001)

但这根本不漂亮。有没有更好的方法来做到这一点？

更新：

我惊讶地发现，如果我将变量设置为 ppm_kind_double 中的常量，问题就不会发生。所以我猜这个“舍入误差”只会发生在我使用的机器的舍入精度位数超过 ppm_kind_double 中定义的数字时。我实际上在集群上运行我的程序（不是这个演示代码），我不知道机器精度。那么也许是那台机器上的四精度导致了问题？

在我将常量设置为双精度后：

before = -0.112_ppm_kind_double
dx = 0.008_ppm_kind_double

【问题讨论】：

你想要-14让我有点失望。我花了一段时间才意识到 -0.112 / 0.008 正好是 -14，但由于舍入误差，结果稍微多一些，然后 CEILING 放大了这个微小的误差。

@chw21，其实我看到了你之前的评论，正在等你发现。 :)

鉴于精度是这里问题的一部分，值得注意的是您的文字常量只是单精度 - 例如-0.112 应该写成-0.112_ppm_kind_double。不过，此更改不会解决根本问题。

除了使常量具有一致的数据类型之外，很难猜出为什么你没有使用 ANINT 而不是 CEILING。

顺便说一句，您的解释是错误的。实际发生的情况是，最初您将变量设置为不准确的单精度，并且使用正确的精度后缀。我强烈建议您阅读 Fortran 中表达式的数值准确性。默认数字是单精度。而你的“机器的四舍五入精度”不存在。

【参考方案1】：

这有点棘手，因为您永远不知道舍入误差从何而来。如果dx 只比0.008 大一点点，那么除法before/dx 可能仍会舍入到相同的值，但现在-13 将是正确的答案。

也就是说，我所见过的最常见的方法是将之前的值轻推到相反的方向。像这样的：

program sign_test
    use iso_fortran_env
    implicit none
    real(kind=real64) :: a, b
    integer(kind=int32) :: c
    a = -0.112
    b = 0.008
    c = my_ceiling(a/b)
    print*, a, b, c
contains
    function my_ceiling(v)
        implicit none
        real(kind=real64), intent(in) :: v
        integer(kind=int32) :: my_ceiling
        my_ceiling = ceiling(v - 1d-6, kind=int32)
    end function my_ceiling
end program sign_test

这不会对绝大多数值产生任何影响，但现在有一些值会被四舍五入超过预期。

【讨论】：

感谢chw21，给初步猜测一个合理的转变也是一个好办法。

【参考方案2】：

请注意，如果您的实数在概念上“精确”到指定的精度，您可能会执行以下操作：

 after=nint(1000*before)/nint(1000*dx) 

这适用于您的示例..您还没有说出您对正值等的期望值，因此您可能需要做一些工作。