随机数生成器和 4 字节与 8 字节整数中的溢出

著名的线性同余随机数生成器也称为最小标准使用公式

x(i+1)=16807*x(i) mod (2^31-1)

我想用 Fortran 来实现这个。

然而，正如《Numerical Recipes》所指出的，直接用默认的Integer类型（32位）实现公式会导致16807*x(i)溢出。

所以本书推荐Schrage的算法是基于m的近似因式分解。该方法仍然可以使用默认的整数类型来实现。

然而，我想知道 fortran 实际上有Integer(8)其范围为的类型-9,223,372,036,854,775,808 to 9,223,372,036,854,775,807这比16807*x(i)可能。

但书上还说了下面这句话

不可能直接实现方程（7.1.2）和（7.1.3） 在高级语言中，由于 a 和 m − 1 的乘积超过了 32 位整数的最大值。

那么为什么我们不能直接使用Integer(8)输入直接执行公式？

是否可以拥有 8 字节整数取决于您的编译器和系统。更糟糕的是传递给的实际值kind获得特定的精度并没有标准化。虽然我知道的大多数 Fortran 编译器都使用字节数（因此 8 就是 64 位），但这并不能保证。

您可以使用selected_int_kind方法得到一种int是有一定范围的。这段代码在我的 64 位计算机上编译并且运行良好：

program ran
    implicit none
    integer, parameter :: i8 = selected_int_kind(R=18)
    integer(kind=i8) :: x
    integer :: i
    x = 100
    do i = 1, 100
        x = my_rand(x)
        write(*, *) x
    end do

    contains
        function my_rand(x)
            implicit none
            integer(kind=i8), intent(in) :: x
            integer(kind=i8) :: my_rand
            my_rand = mod(16807_i8 * x, 2_i8**31 - 1)
        end function my_rand
end program ran

以下@VladimirF 评论的更新和解释

现代 Fortran 提供了一个名为iso_fortran_env提供引用标准变量类型的常量。在你的情况下，人们会使用这个：

program ran
    use, intrinsic :: iso_fortran_env, only: int64
    implicit none
    integer(kind=int64) :: x

然后如上所述。这段代码比旧的更容易阅读selected_int_kind。 （为什么R一定要回到18岁吗？）