Universidad de Costa Rica Centro de Investigación en Matemática Pura y Aplicada

Park-Miller-Carta Random Number Generator in Fortran

Abstract   

Motivation   

When Park & Miller developed their "minimal standard" for a linear congruential [LCG] pseudo random number generator [PRNG] they succeeded in providing a "good enough" LCG for "the rest of us" [PM-1998]. In their paper, they provided a Pascal implementation that can be easily translated into C or C++, as the basic formula for the LCG is simple: (16807 mod 2^32 -1). Later on, Carta suggested and improvement of the algorithm that uses no division [Carta-1990]. A good explanation of Carta's improvement can be found in [Whittle-2013], where a Fortran implementation was promised. That promise gets fulfilled in here.

#include <stdint.h>
int32_t next_16807_64b( int32_t seed ) {
    enum { two_31_1 = 2147483647 };
    int64_t next = seed;
    next = ( ( next * (int64_t)(16807) ) % two_31_1 );
    return (int32_t)(next);
}

Figure 1 is a straightforward implementation of the "minimal standard". It has two defects: uses 64 bit arithmetic and division. This LCG will generate every number in range "[1,...,2^31-2]" exactly once: this is a full period generator.

Schrage's improvement   

In [Schrage-1987] a simple procedure is provided to calculate the result of a multiplication using modulo arithmetic. This trick is based on the Quotient/Remainder theorem that states that for any integer "m" and "a>0" there exist unique integers "q" (quotient) and "r" (remainder) such that "m = a*q+r" and "0≤r<m".

This theorem was applied by Schrage to avoid overflow when calculating the modulo result of an integer multiplication [Bathsheba-2014] . Suppose you want to calculate the multiplication "a*z mod m". First, decompose the modulus "m>0" using the multiplier "a": "m = a*q+r". This means that "q=[m/a]" and "r = m mod a", where [] is the integer part. It does not always happen, but suppose that "r<q"and "0<z<m-1" then both "a*(z mod q)" and "r*[z/q]" will be in range "[0,...,m-1]". Furthermore, "a*z mod m = a*(z mod q) - r*[z/q]". If this difference is negative just add "m".

function Random32 : real;
(* Integer Version 2 *)
const
    a = 16807;
    m = 2147483647;
    q = 127773; (* m div a *)
    r = 2836;   (* m mod a *)
var
    lo, hi, test : integer;
begin
    hi := seed div q;
    lo := seed mod q;
    test := a * lo - r * hi;
    if test > 0 then
        seed := test
    else
        seed := test + m;
    Random32 := seed / m
end;

The implementation in Figure 2 is the original implementation provided by Park & Miller in their paper, where they use Schrage's trick to calculate "a*z mod m". No overflow will occur and 32 bit arithmetic is used, but division is still used in the algorithm.

Carta's implementation   

To avoid using divisions, Carta's implementation of the "minimal standard" uses only multiplications [Carta-1990]. His explanation is a little difficult to follow, but in a latter paper Whittle provides a simpler explanation [Whittle-2013].

/* Park-Miller "minimal standard" 31 bit
 * pseudo-random number generator, implemented
 * with Davis G. Carta's optimisation: with
 * 32 bit math and without division.
 */

long unsigned int rand31_next()
{
    long unsigned int lo, hi;

    lo = 16807 * (seed & 0xFFFF);
    hi = 16807 * (seed >> 16);

    lo += (hi & 0x7FFF) << 16;
    lo += hi >> 15;

    if (lo > 0x7FFFFFFF) lo -= 0x7FFFFFFF;

    return ( seed = (long) lo );
}

The implementation in Figure 3 is taken from Whittle's contribution. It no longer uses division and relies on 32 bit unsigned arithmetic. To test that it does produce the correct result, its output can be compared with the original algorithm in Figure 2; on current machines, this takes little time (less than a few minutes). 2's complement arithmetic is assumed in this implementation (a ternary computer might not produce the correct result).

Fortran issues   

When you have years of C++ programming, it becomes hard to produce Fortran implementations. However, it is easy enough to translate Carta's C implementation into Fortran, but small difficulties arise. When testing whether variable "lo" is negative, the C language test uses a comparison to "0x7FFFFFFF" ( 2^31-1 ). As Fortran integers are signed the unsigned integer test in C cannot be directly transliterated to Fortran: another test that achieves the same result must be used.

      subroutine pmc_next( SEED )
         integer(4) SEED
         integer(4) lo, hi

         lo = 16807 * IAND( SEED, 65535_4 )   ! Z'FFFF' = 65535
         hi = 16807 * ISHFT(SEED,-16)
         lo = lo + ISHFT(IAND(hi,32767_4),16) ! Z'7FFF' = 32767
         lo = lo + ISHFT(hi,-15)

         if ( lo .LT. 0 ) then    ! 2147483647==Z'7FFFFFFF'
            lo = lo-2147483647_4
         endif

         SEED    = lo
         return
      end subroutine

When 2's complement is used in computers, the first bit of an integer number would be '1' when the stored value is negative. Hence, instead of asking whether the "0x7FFFFFFF" barrier was crossed, the signed equivalent is to test for a negative value. This is apparent in the Fortran implementation of Figure 4. The underscored suffix digit in Fortran denotes the constant byte size: "2147483647_4" is a 4 byte integer value (32 bits) and "2147483647_8" is an 8 byte integer value (64 bits).

      Program test_rand31pmc
c     ----------------------

      integer(4) SEED
      integer(4) next_fortran,  next_cpp

      write (*,'(A)') 'BEGIN'
      do SEED = 1 , 2147483647_4-1
         next_fortran = SEED;
         call pmc_next( next_fortran );

         next_cpp = SEED;
         call pm_next( next_cpp );

         if ( next_fortran .NE. next_cpp ) then
            write (*,'(I10,A,I0,A,I0)')
     c         SEED,' -> ',next_cpp,'.NE.',next_fortran 
               call sleep(1*0)
         endif
      enddo

      write (*,'(A)') 'END'

      if (.true.) then
         write (*,'(I0,A,I10,A)')
     c      sizeof(2147483647_4),'->int(4)',
     c      sizeof(2147483647_8),'->int(8)'
      endif

      end ! test_rand31pmc

Testing the Fortran implementation for correctness is done with a Fortran main program linked with the C language version of the Park & Miller algorithm. A simple loop through every possible value test that the Fortran implementation always produces the correct value, as shown in Figure 5.

Some other issues had to be tackled to mix Fortran and C, but the problems rise when using an Integrated Development Environment [IDE]. In this case Code::Blocks version 16.01 was used, but it is a little stiff for mixing modules from different languages. In the Windows platform, the Code::Blocks Fortran project file (with extension .cbp) required to use a direct path to invoke the C compiler:

<Option
  compiler="gfortran"
  use="1"
  buildCommand='&quot;%ProgramFiles%\CodeBlocks\MinGW\bin\gcc.exe&quot; -g -c $file -o $object'
/>

Fewer problems arise when invoking the Fortran and C compilers from the command line, as Unix/Linux people usually do. A good guide with more information on mixed language programming is [Ippolito-2016].

Conclusions   

The "minimal standard" LCG has many qualities. Implementation that use Carta's speed up help choosing it for many applications where "cryptographically secure" PRNG are not required.

Acknowledgments   

      Alejandro Di Mare made valuable suggestions that helped improve earlier versions of this work. The Centro de Investigación en Matemática Pura y Aplicada [CIMPA] at the Universidad de Costa Rica provided resources for this research. The spelling and grammar were checked using the http://spellcheckplus.com/ tool.

References   

Index   

[1]	Motivation
[2]	Schrage's improvement
[3]	Carta's implementation
[4]	Fortran issues
[5]	Conclusions
[6]	Acknowledgements
	References
	Index
	About the author
	About this document
	Begin Index End

About the author   

Adolfo Di Mare <adolfo@di-mare.com>

About this document   

Referencia:	Di Mare, Adolfo: Park-Miller-Carta Random Number Generator in Fortran: Technical Report 2016-02-ADH, Centro de Investigación en Matemática Pura y Aplicada [CIMPA], Universidad de Costa Rica, 2016.
Internet:	http://www.di-mare.com/adolfo/p/rand31pmc.htm       Google™ Translate http://www.di-mare.com/adolfo/p/src/rand31pmc.zip
Author:	Adolfo Di Mare <adolfo@di-mare.com>
Contact:	Apdo 4249-1000, San José Costa Rica Tel: (506) 2511-6606       Fax: (506) 2511-4918
Revision:	CIMPA, September 2016
Visitantors: