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?

`Integer(8)` does not mean 8-byte integer, but integer of kind 8. It is not the same http://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 — Vladimir F Героям слава, Jul 21 '16 at 10:01
Regarding Numerical Recipes see http://www.uwyo.edu/buerkle/misc/wnotnr.html . Simply don't believe them for this. It is an old book and wasn't that good even when it was written. — Vladimir F Героям слава, Jul 21 '16 at 10:05

chw21 · Accepted Answer · 2022-11-21T09:46:32.583

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?)

Thank you chw21. So the saying of the book is wrong, the formula can be directly implemented in fortran, right? — user15964, Jul 21 '16 at 09:16
@user15964 Yes, you can. To select 8 bytes directly you can use `kind=int64` where `int64` is defined in intrinsic module `iso_fortran_env`. — Vladimir F Героям слава, Jul 21 '16 at 10:02

Anthony Scemama · Answer 2 · 2016-07-21T12:13:52.167

1

Yes. The simplest thing is to append _8 to the integer constants to make them 8 bytes. I know it is "old style" Fortran but is is portable and unambiguous.

By the way, when you write:

16807*x mod (2^31-1)

this is equivalent to take the result of 16807*x and use an and with a 32-bit mask where all the bits are set to one except the sign bit. The efficient way to write it by avoiding the expensive mod functions is:

iand(16807_8*x, Z'7FFFFFFF')

Update after comment :

or

iand(16807_8*x, 2147483647_8)

if your super modern compiler does not have backwards compatibility.

edited Jul 21 '16 at 12:13

answered Jul 21 '16 at 09:42

Anthony Scemama

1,563
12
19

It is not portable. `_8` is not accepted by certain compilers which have kinds `1,2,3`. `_8` DOES NOT mean 8 bytes, it means `kind=8`. http://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 And it is not old style (it only works in Fortran 90 and newer), it is just bad style. – Vladimir F Героям слава Jul 21 '16 at 09:59
I am also strongly convinced there are problems with your `iand(16807_8*x, Z'7FFFFFFF')`. The usage of BOZ constants were changed in Fortran 2008 and are more general now but still I am quite sure that the BOZ constant in your example is quadruple precision (in my compiler) and mixing it this way is not allowed. – Vladimir F Героям слава Jul 21 '16 at 10:25
OK. Thanks for correcting. Fortran seems to deviate so much from Fortran90! It reminds me of the Python2 vs Python3 conflict... – Anthony Scemama Jul 21 '16 at 12:15
No, Fortran 90 did not allow a BOZ constant in your expression at all. It is completely invalid there. Fortran 2008 just allowed them to appear in more contexts. – Vladimir F Героям слава Jul 21 '16 at 12:20
2

As of Fortran 2008, BOZ literal constant in a reference to `iand` is converted to a value that has the same kind as the other argument (both arguments may not be BOZ literals). That reference looks ok to me. BOZ literals do not have a type or kind in their own right in the standard language - that is always inferred from the contexts in which they can be used. – IanH Jul 21 '16 at 12:33
1

Please don't use hard-coded kind numbers such as _8. Not all compilers use 8 as the kind for this integer size. Define a PARAMETER constant, such as I64, with the value SELECTED_INT_KIND(15). Then use _I64 as the kind specifier. – Steve Lionel Jul 21 '16 at 14:16
@IanH I know that they do not have a kind per se. Could you show me the standard reference? Perhaps a bug should be reported to gfortran then. – Vladimir F Героям слава Jul 21 '16 at 21:40
@VladimirF F2008 C4102, plus the various descriptions of what happens where BOZ literals are permitted. – IanH Jul 22 '16 at 23:03

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

2 Answers2

Linked

Related