C:/Musimathics_local/Musimat/MusimatChapter9/C090703.cpp

Go to the documentation of this file.

#include "MusimatChapter9.h"
MusimatChapter9Section(C090703) {
Print("*** 9.7.3 Linear Congruential Method ***");
/*****************************************************************************

9.7.3 Linear Congruential Method

Equation (9.1) shows the linear congruential method for generating random numbers, 
introduced by D. H. Lehmer in 1948 (Knuth 1973, vol. 2):

        x[n+1] = ((ax[n]+b))c,  n>=0.                                                                           (9.1)

The notation ((x))n means “x is reduced modulo n.” The result is the remainder 
after integer division of x by n (see appendix A).
Equation (9.1) is a recurrence relation because the result of the previous step 
(x[n]) is used to calculate a subsequent step (x[n+1]). It is linear because the 
ax + b part of the equation describes a straight line that intersects the y-axis 
at offset b with slope a. Congruence is a condition of equivalence between two 
integers modulo some other integer, and refers here simply to the fact that 
modulo arithmetic is being used.  

For successive computations of x, the output will grow until it reaches the value c. 
When c is exceeded, the new value of x is effectively reset to x – c by the modulus operation. 
A new slope will grow from this point, and this process repeats endlessly.
The result can be quite predictable depending upon the values of a, b, c, and x[0],
the initial value of x.

For instance, if a = b = x[0] = 1, and c = infinity, an ascending straight line at a 45 degree 
slope is produced. However, for other values, the numbers generated can appear random.
In practice, the modulus c should be as large as possible in order to produce long 
random sequences. On a computer, the ultimate limit of c is the arithmetic precision 
of that machine. For example, if the computer uses 16-bit arithmetic, random numbers 
generated by this method can have at most a period of 2^16 = 65,536 values before 
the pattern repeats.

The quality of randomness within a period varies depending on the values chosen 
for a, x, and b. Much heavy-duty mathematics has been expended choosing good 
values (Knuth 1973, vol. 2). For 32-bit arithmetic, Park and Miller (1988) recommend 
a = 16,807, b = 0, and c = 2,147,483,647.
The linear congruential method is appealing because once a good set of the parameters is 
found, it is very easy to program on a computer. The LCRandom() method returns a random 
number by the linear congruential method each time it is called:
*****************************************************************************/
        para1(); // Step into this function to continue.
}

//Constants from Park and Miller
Const Integer a = 16807;                        // a, b, c and x are global constant values
Const Integer b = 0;
Const Integer c = 2147483647;
Integer x = 1;                                          // x stores the value produced by
                                                                        // LCRandom between invocations
Integer LCRandom(){
        x = Mod(a * x + b, c);                  // update x based on its previous value
        Integer r = x;                                  // x may be positive or negative
        If (r < 0)                                              // force the result to be positive
                r = -r;
        Return(r);
}


Static Void para1() {
/*****************************************************************************
The parameters a, b, and c are constant (time-invariant) system parameters. 
Parameter x is initialized in this example to 1, but it can be initialized 
to any other integer. The value of a * x + b is calculated, the remainder 
is found modulo c, and the result is reassigned to x.

While the value of x is less than c, x grows linearly. When the expression 
a * x + b eventually produces a value beyond the range of c, then x is reduced 
modulo c. The random effect of this method comes from the surprisingly unpredictable 
sequence of remainders generated by the modulus operation, depending upon careful 
choice of parameters.

The calculation of x ranges over all possible positive and negative integers 
smaller than the value of ±c. But it is generally preferable to constrain its 
choices to a range. To make this conversion easier, we force the result to be a 
positive integer.

Below is a sample invocation of LCRandom().
*****************************************************************************/
        Print("*** Ten invocations of LCRandom() ***");
        IntegerList x;
        
        For (Integer i = 0; i < 10; i++ )
                x[i] = LCRandom();

        Print( x );
}}

/* $Revision: 1.4 $ $Date: 2006/09/12 17:37:24 $ $Author: dgl $ $Name:  $ $Id: _c090703_8cpp-source.html,v 1.4 2006/09/12 17:37:24 dgl Exp $ */
// The Musimat Tutorial © 2006 Gareth Loy
// Derived from Chapter 9 and Appendix B of "Musimathics Vol. 1" © 2006 Gareth Loy 
// and published exclusively by The MIT Press.
// This program is released WITHOUT ANY WARRANTY; without even the implied 
// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
// For information on usage and redistribution, and for a DISCLAIMER OF ALL
// WARRANTIES, see the file, "LICENSE.txt," in this distribution.
// "Musimathics" is available here:     http://mitpress.mit.edu/catalog/item/default.asp?ttype=2&tid=10916
// Gareth Loy's Musimathics website:    http://www.musimathics.com/
// The Musimat website:                 http://www.musimat.com/
// This program is released under the terms of the GNU General Public License
// available here:                      http://www.gnu.org/licenses/gpl.txt

---