00001 #include "MusimatChapter9.h" 00002 MusimatChapter9Section(C090703) { 00003 Print("*** 9.7.3 Linear Congruential Method ***"); 00004 /***************************************************************************** 00005 00006 9.7.3 Linear Congruential Method 00007 00008 Equation (9.1) shows the linear congruential method for generating random numbers, 00009 introduced by D. H. Lehmer in 1948 (Knuth 1973, vol. 2): 00010 00011 x[n+1] = ((ax[n]+b))c, n>=0. (9.1) 00012 00013 The notation ((x))n means "x is reduced modulo n." The result is the remainder 00014 after integer division of x by n (see appendix A). 00015 Equation (9.1) is a recurrence relation because the result of the previous step 00016 (x[n]) is used to calculate a subsequent step (x[n+1]). It is linear because the 00017 ax + b part of the equation describes a straight line that intersects the y-axis 00018 at offset b with slope a. Congruence is a condition of equivalence between two 00019 integers modulo some other integer, and refers here simply to the fact that 00020 modulo arithmetic is being used. 00021 00022 For successive computations of x, the output will grow until it reaches the value c. 00023 When c is exceeded, the new value of x is effectively reset to x - c by the modulus operation. 00024 A new slope will grow from this point, and this process repeats endlessly. 00025 The result can be quite predictable depending upon the values of a, b, c, and x[0], 00026 the initial value of x. 00027 00028 For instance, if a = b = x[0] = 1, and c = infinity, an ascending straight line at a 45 degree 00029 slope is produced. However, for other values, the numbers generated can appear random. 00030 In practice, the modulus c should be as large as possible in order to produce long 00031 random sequences. On a computer, the ultimate limit of c is the arithmetic precision 00032 of that machine. For example, if the computer uses 16-bit arithmetic, random numbers 00033 generated by this method can have at most a period of 2^16 = 65,536 values before 00034 the pattern repeats. 00035 00036 The quality of randomness within a period varies depending on the values chosen 00037 for a, x, and b. Much heavy-duty mathematics has been expended choosing good 00038 values (Knuth 1973, vol. 2). For 32-bit arithmetic, Park and Miller (1988) recommend 00039 a = 16,807, b = 0, and c = 2,147,483,647. 00040 The linear congruential method is appealing because once a good set of the parameters is 00041 found, it is very easy to program on a computer. The lcRandom() method shown below returns a random 00042 number by the linear congruential method each time it is called. lcRandom() is a 00043 copy of the Musimat LCRandom() routine in Random.cpp reproduced here for 00044 didactic reasons (to keep things simple). 00045 *****************************************************************************/ 00046 para1(); // Step into this function to continue. 00047 } 00048 00049 //Constants from Park and Miller 00050 Const Integer lc_a = 16807; // a, b, c and x are global constant values 00051 Const Integer lc_b = 0; 00052 Const Integer lc_c = 2147483647; 00053 Integer lc_x = 1; // x stores the value produced by LCRandom between invocations 00054 00055 Integer lcRandom(){ 00056 lc_x = Mod(lc_a * lc_x + lc_b, lc_c); // update x based on its previous value 00057 Integer r = lc_x; // x may be positive or negative 00058 If (r < 0) // force the result to be positive 00059 r = -r; 00060 Return(r); 00061 } 00062 00063 00064 Static Void para1() { 00065 /***************************************************************************** 00066 The parameters a, b, and c are constant (time-invariant) system parameters. 00067 Parameter x is initialized in this example to 1, but it can be initialized 00068 to any other integer. The value of a * x + b is calculated, the remainder 00069 is found modulo c, and the result is reassigned to x. 00070 00071 While the value of x is less than c, x grows linearly. When the expression 00072 a * x + b eventually produces a value beyond the range of c, then x is reduced 00073 modulo c. The random effect of this method comes from the surprisingly unpredictable 00074 sequence of remainders generated by the modulus operation, depending upon careful 00075 choice of parameters. 00076 00077 The calculation of x ranges over all possible positive and negative integers 00078 smaller than the value of ±c. But it is generally preferable to constrain its 00079 choices to a range. To make this conversion easier, we force the result to be a 00080 positive integer. 00081 00082 Below is a sample invocation of LCRandom(). 00083 *****************************************************************************/ 00084 Print("*** Ten invocations of LCRandom() ***"); 00085 IntegerList x; 00086 00087 For (Integer i = 0; i < 10; i++ ) { 00088 x[i] = lcRandom(); 00089 } 00090 00091 Print( x ); 00092 } 00093 00095 /* $Revision: 1.3 $ $Date: 2006/09/05 08:02:45 $ $Author: dgl $ $Name: $ $Id: C090703.cpp,v 1.3 2006/09/05 08:02:45 dgl Exp $ */ 00096 // The Musimat Tutorial © 2006 Gareth Loy 00097 // Derived from Chapter 9 and Appendix B of "Musimathics Vol. 1" © 2006 Gareth Loy 00098 // and published exclusively by The MIT Press. 00099 // This program is released WITHOUT ANY WARRANTY; without even the implied 00100 // warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 00101 // For information on usage and redistribution, and for a DISCLAIMER OF ALL 00102 // WARRANTIES, see the file, "LICENSE.txt," in this distribution. 00103 // "Musimathics" is available here: http://mitpress.mit.edu/catalog/item/default.asp?ttype=2&tid=10916 00104 // Gareth Loy's Musimathics website: http://www.musimathics.com/ 00105 // The Musimat website: http://www.musimat.com/ 00106 // This program is released under the terms of the GNU General Public License 00107 // available here: http://www.gnu.org/licenses/gpl.txt
1.7.2