Overflow in a random number generator and 4-byte vs. 8-byte integers

Question

The famous linear congruential random number generator also known as minimal standard use formula

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

I want to implement this using Fortran.

However, as pointed out by "Numerical Recipes", directly implement the formula with default Integer type (32bit) will cause 16807*x(i) to overflow.

So the book recommend Schrage’s algorithm is based on an approximate factorization of m. This method can still implemented with default integer type.

However, I am wondering fortran actually has Integer(8) type whose range is -9,223,372,036,854,775,808 to 9,223,372,036,854,775,807 which is much bigger than 16807*x(i) could be.

but the book even said the following sentence

It is not possible to implement equations (7.1.2) and (7.1.3) directly in a high-level language, since the product of a and m − 1 exceeds the maximum value for a 32-bit integer.

So why can't we just use Integer(8) type to implement the formula directly?

Answer 1

Whether or not you can have 8-byte integers depends on your compiler and your system. What's worse is that the actual value to pass to kind to get a specific precision is not standardized. While most Fortran compilers I know use the number of bytes (so 8 would be 64 bit), this is not guaranteed.

You can use the selected_int_kindmethod to get a kind of int that has a certain range. This code compiles on my 64 bit computer and works fine:

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

Update and explanation of @VladimirF's comment below

Modern Fortran delivers an intrinsic module called iso_fortran_env that supplies constants that reference the standard variable types. In your case, one would use this:

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

and then as above. This code is easier to read than the old selected_int_kind. (Why did R have to be 18 again?)