Math 455.1  •  29 March 2009

Linear congruential generator lrand48

Copyright © 2009 by Murray Eisenberg.  All rights reserved.

The algorithm

This notebook illustrates how congruences may be used to generate pseudorandom integers. The algorithm is used by the linear congruential pseudorandom number generator function lrand48() that forms part of the standard C library. Here is the description of the function, as extracted by issuing in Linux the command man lrand48.

All the routines work by generating  a  sequence  of  48-bit
     integer, Xi , according to the linear congruential formula

          Xn+1 = (aXn + c) mod m, where n >= 0.

     The parameter m=2**48, hence 48-bit  integer  arithmetic  is
     performed. Unless  lcong48() is called, a and c are given by:

           a = 0x5DEECE66D
           c = 0xB

In that description, Xi, Xn+1, and Xn really mean subscripted quantities (Subscript[X, i]), (Subscript[X, n+1]), and (Subscript[X, n]). And the forms 0x... given for the numbers a and c are their hexadecimal (base 16) representations.

Implement the algorithm in Mathematica

Set up the coefficients:

aOctal=273673163155;
a=FromDigits[IntegerDigits[aOctal],8]

cOctal=13;
c=FromDigits[IntegerDigits[cOctal],8]

The modulus—and the number of distinct pseudorandom integers generated by this algorithm—is:

m=(2^48)

The sequence of pseudorandom integers will start with—be "seeded by"—the sum of the entries in the current date.

DateList[]

The next cell defines a function x whose values are successive pseudorandom integers starting with the random seed x[0].

Clear[x]
x[n_]:=x[n]=Mod[a x[n-1]+c,m]
x[0]=Plus@@Round@DateList[];

Note.  The expression x[n_] := x[n] = in the preceding definition caches the values as they are obtained.

Use the Mathematica implementation

Table[x[n],{n,0,20}]//TableForm