Fortran 90 ran_init 整数模型假设中的数值方法

Question

数值食谱的 ran_init 子例程包含以下行：

INTEGER(K4B) :: new,j,hgt
...                                                                                            
hgt=hg 
...
if (hgt+1 /= hgng)    call nrerror('ran_init: arith assump 3 fails')

其中 K4B、hgng 和 hg 在模块中通过以下方式全局声明：

INTEGER, PARAMETER :: K4B=selected_int_kind(9) 
INTEGER(K4B), PARAMETER :: hg=huge(1_K4B), hgm=-hg, hgng=hgm-1

问题是在一台特定的计算机上（但在其他计算机上没有）出现错误 ran_init: arith assump 3 fails。我从文档中得到的关于此错误的唯一信息是：

Bit of dirty laundry here! We are testing whether the most positive integer hg wraps around to the most negative integer hgng when 1 is added to it. We can’t just write hg+1 , since some compilers will evaluate this at compile time and return an overflow error message. If your compiler sees through the charade of the temporary variable hgt , you’ll have to find another way to trick it!

我该如何欺骗它？

Answer 1

崩溃的直接原因：

编译器在优化步骤中可以很容易地证明 hgt+1 和 hgng 不能相同，因此它将整个条件优化为 .false.。两个正整数相加不能得到一个负整数。

您可以通过使用较小级别的优化或使用一些临时变量的伪装来避免它，这些伪装比他们在旧的数字食谱中所做的更复杂。您可以尝试使用 -fwrapv 或 -fno-strict-overflow 标志。但它非常不确定，只能在具有特定编译器标志的给定版本的编译器中工作。

最简单的 "fix" 可能只是删除有问题的检查。它很可能在许多情况下都有效（而在其他情况下则非常失败）。

解释及解决办法：

Numerical Recipes 所做的是违反 Fortran 标准的。他们假设如果将最大整数加 1，就会得到最小整数。

这不允许用于 Fortran 整数，也不允许用于 C 中的有符号整数。它只允许用于 C 中的无符号整数，而 Fortran 没有。

因此编译器可以安全地假设将两个正整数相加永远不会得到负整数。有关这些优化的详尽讨论，请参阅 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=30475。

一个选择是重写C中可能导致溢出的操作，并使用转换为无符号整数。我在我使用的随机数生成器中使用它（基于http://www.cmiss.org/openCMISS/wiki/RandomNumberGenerationWithOpenMP）：

#include <stdint.h>
#include <memory.h>

int32_t sum_and_overflow (int32_t a, int32_t b)
{
  uint32_t au, bu, su;
  int32_t s;

  memcpy(&au, &a, sizeof(a));
  memcpy(&bu, &b, sizeof(b));
  su = au + bu;
  memcpy(&s, &su, sizeof(s));
  return s;
}

使用 Fortran 界面

interface
  function sum_and_overflow(a, b) result(res) bind(C, name="sum_and_overflow")
    use, intrinsic :: iso_c_binding
    integer(c_int32_t) :: res
    integer(c_int32_t), value :: a, b
  end function
end interface

而不是

c = a + b

我打电话

c =  sum_and_overflow(a, b)

所以在你的情况下

if (sum_and_overflow(hgt,1) /= hgng) ...

但不仅如此，您还必须至少在一个地方找到在生成器中使用此假设的地方。在我使用的生成器中只有一根这样的线。

还有许多其他 hack 可能会导致暂时成功，但稍后会因其他一些编译器选项而失败。例如，GCC 中的未定义行为清理不喜欢 Fortran 和 C 中的有符号整数溢出，如果您这样做，将终止您的程序。

这就是为什么我尝试展示一个更复杂的解决方案，但不是围绕标准工作，而是遵循它。

Fortran 90 ran_init 整数模型假设中的数值方法

Numerical Recipes in fortran 90 ran_init integer model assumption

fortran

integer

numerical-recipes

崩溃的直接原因：

解释及解决办法：