[UCR]  
[/\]
Universidad de Costa Rica
Centro de Investigación en
Matemática Pura y Aplicada
[<=] [home] [<>] [\/] [=>]
Google Translate

Park-Miller-Carta Random Number Generator in Fortran

Adolfo Di Mare




Abstract [<>] [\/] [/\]

The Park and Miller "minimum standard" pseudo random number generator is widely used because of its simple implementation. The improvement suggested by Carta yields a faster linear congruential generator that uses no division. The Fortran implementation provided here complements the already available C and C++ versions. El "estándar mínimo" para generar números seudo aleatorios desarrollado por de Park y Miller es ampliamente utilizado debido a su sencilla implementación. La mejora sugerida por Carta produce un generador congruencial lineal más rápido que no utiliza división. La implementación Fortran aquí proporcionada complementa las ya disponibles en C y C++.

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

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;
Figure 2

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 );
}
Figure 3

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
Figure 4

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
Figure 5

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 [<>] [\/] [/\]

[Bathsheba-2014] Bathsheba: Modular multiplication of large numbers in c++, 2014.
      http://stackoverflow.com/questions/20971888/#20972369
[Carta-1990] Carta, David G.: Two Fast Implementations of the "Minimal Standard" Random Number Generator, 1990.
      http://www.firstpr.com.au/dsp/rand31/p87-carta.pdf
[Ippolito-2016] Ippolito, Greg: Using C/C++ and Fortran together, 2001-2016.
      http://www.yolinux.com/TUTORIALS/LinuxTutorialMixingFortranAndC.html
[PM-1998] Park, Stephen K. & Miller, Keith W.: Random Number Geuerators: Goodones Are Hard To Find, Communications of the ACM, Oct 1988, Vol 31 Number 10 1192-1201.
      http://www.firstpr.com.au/dsp/rand31/p1192-park.pdf
[Schrage-1987] Bratley, Paul & Fox, Bennet L. & Schrage, Linus E.: A Guide to Simulation 2ed, ISBN-10: 0387964673, ISBN-13: 978-0387964676, Springer 1987.
[Whittle-2013] Whittle, Robin: Park-Miller-Carta Pseudo-Random Number Generator, 2005-2013.
      http://www.firstpr.com.au/dsp/rand31

           

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: Investigador costarricense en el Centro de Investigación en Matemática Pura y Aplicada [CIMPA] y en la Escuela de Ciencias de la Computación e Informática [ECCI] de la Universidad de Costa Rica [UCR], en donde ostenta el rango de Profesor Catedrático. Trabaja en las tecnologías de Programación e Internet. También es Catedrático de la Universidad Autónoma de Centro América [UACA]. Obtuvo la Licenciatura en la Universidad de Costa Rica, la Maestría en Ciencias en la Universidad de California, Los Angeles [UCLA], y el Doctorado (Ph.D.) en la Universidad Autónoma de Centro América.
Adolfo Di Mare: Costarrican Researcher at the Centro de Investigación en Matemática Pura y Aplicada [CIMPA] and the Escuela de Ciencias de la Computación e Informática [ECCI], Universidad de Costa Rica [UCR], where he is full professor and works on Internet and programming technologies. He is Cathedraticum at the Universidad Autónoma de Centro América [UACA]. Obtained the Licenciatura at UCR, and the Master of Science in Computer Science from the University of California, Los Angeles [UCLA], and the Ph.D. at the Universidad Autónoma de Centro América.
[mailto]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:

 

Copyright © 2016 Adolfo Di Mare
Derechos de autor reservados © 2016 Adolfo Di Mare <adolfo@di-mare.com>
[home] [<>] [/\]