Talk:Linear congruential generator
From Wikipedia, the free encyclopedia
Contents[hide] |
[edit] C code example of LCG for 32-bit CPU
C code example of LCG for 32-bit (or more) CPU:
const unsigned M = 0xffffffff - 8; /* 2^32 - 9 = 4294967287 = 13 * 71^2 * 65539 */
const unsigned A = 13 * 71 * 65539 + 1;
const unsigned B = 17; /* or any other relatively prime to M number */
int main() {
unsigned c = 0;
unsigned v = 1;
do {
v = ((int64_t)A * v + B) % M;
c++;
if (!(c & 0x3ffffff))
printf("%10u: %08x\n", c, v);
} while (v != 1);
printf("%u (0x%x) iterations: loop detected (value %08x repeats)\n", c, c, v);
return 0;
}
Output:
67108864: 7a15a240
134217728: 63f193d7
...
4227858432: c3361633
4294967287 (0xfffffff7) iterations: loop detected (value 00000001 repeats)
Is it worth adding to the article?
[edit] LCG examples
The Mersenne twister both runs faster than and generates higher-quality deviates than almost any LCG...
While not crypto-strong enough, LCG is fast to set up and calculate, and needs just one word of storage. For many applications this is good enough (why should I bother with creating array of 624 words for Mersenne twister if I only need a 8-byte salt for my new Unix password?). Yet the article contains only one actual example of generator (A = 1664525, B = 1013904223, M = 2^32), which is particularly bad because of that alternating even-odd sequence.
http://www.gnu.org/software/gsl/manual/html_node/Other-random-number-generators.html mentions some other LCGs. Is it worth adding a few of those as "good" examples too, so that people do not pick parameters which are particularly bad?
x_{n+1} = (16807 * x_n) mod (2^31 - 1) - period is 2^31-1 (right?)
x_{n+1} = (48271 * x_n) mod (2^31 - 1) - period is ???
x_{n+1} = (62089911 * x_n) mod (2^31 - 1) - period is ???
- The above three generator have period 231 − 2, the maximum period for generators of this type. Maxal (talk) 23:09, 22 February 2008 (UTC)
x_{n+1} = (40692 * x_n) mod (2^31 - 249) - period is ???
- This one has period 231 − 250, which is also the maximum one for Park-Miller RNG modulo prime 231 − 249. Maxal (talk) 23:09, 22 February 2008 (UTC)
x_n = (271828183 * x_{n-1} + 314159269 * x_{n-2}) mod (2^31 - 1)
which isn't a LCG. Is it "better" or what?
- It is a linear congruent sequence of the second order. It may be better is a sense that its period may be equal m2 not just m as for linear congruent sequence of the first order. Maxal (talk) 22:58, 22 February 2008 (UTC)
Gotta mention the classic Speccy one; x_{n+1} = (75 * (x_n + 1) - 1) mod (2^16 + 1) - period is 2^16.
- This type of LCGs (i.e., with c=0) is a different beast. I've separated it in Park-Miller RNG article. Speccy's one will go there. Maxal (talk) 22:58, 22 February 2008 (UTC)
[edit] Grammar issues
Who ever wrote this article is obviously good at mathematics but has very poor grasp of using plain language. The grammar is extremely bad making the article baffling to anybody but the already knowledgable.
- Yes, it's a lovely article kids, but for those of us on planet Earth a little more explanation might be a good idea. Or not, as you like. Greglocock 11:15, 21 September 2006 (UTC)
[edit] disproof
The following disproof is not valid:
"Choose M=11, A=5, B=1, V=1 All conditions are satisfied:-
1) 1 and 11 are relatively prime. 2) 11 is prime and has no prime factors. 3) A-1 = 4 , 4) M > max(5, 1, 1) 5) A > 0 , B > 0
)
Yet this yields the following sequence:- 6 9 2 0 1 6 9 2 0 1 6 "
It fails because when 11 is reduced to prime factors, the only prime factor is 11. A-1=4 is not divisible by 11, so this set of inputs fails the 2nd rule. See the article on prime factor.
[edit] Rules for the maximum period
Those rules can be found in the following literature: A. M. Law and W. David Kelton, Simulation Modeling and Analysis, McGraw-Hill, New York, 1991
T. Hull and A. Dobell, Random Number Generators, SIAM Rev., 4: 230-252, 1962 ( http://locus.siam.org/SIREV/volume-04/art_1004061.html )
The text found in Hull and Dobell states that "The sequence defined by the congruence relation (1) has full period m, provided that: (i) c is relatively prime to m; (ii) a = 1 (mod p) if p is a prime factor of m; (iii) a = 1 (mod 4) if 4 is a factor of m. ", indicating that the conditions are sufficient but not necessary. (note that the formula they use is labeled as "ax + c (mod m)" ).
It seems that these rules do not work when M is prime, given what was said above in this discussion page, presumeably because of the 2nd rule. Given that, the conditions being sufficient but not necessary seems to be a popular mis-citation. I was unable to find clarifications of the rules (they are not as precice as they should be in the original paper). If there is a flaw with these rules, or if the original proof has been disproven, then many students have been taught these rules in error. Many citations in lecture notes can be found of them by searching for the literal string in www.google.com:
"If q is a prime number that divides m, then q divides a-1"
Even if the only mis-citation is that the rules are necessary and sufficient, many students are still being taught incorrectly.
[edit] Case B=0
Something like the following is needed:
The maximal period of V_{j+1} = A×V_j + B mod M, is M-1 if B = 0 (0 is fix point), otherwise it is M.
Case 1. If B = 0 ⇒ M must be prime for a full period (otherwise starting with a proper divisor V_0 of M, all V_j values are divisors).
V_{j+1} = A×V_j mod P, with P prime, is of full period (cycle length = P–1), if and only if P divides A^i − 1 for i = P−1, but for no smaller i.
Case 2. (B != 0) V_{j+1} = ( A×V_j + B ) mod M is full period
if B and M are relatively prime, and if a prime p divides M, then p divides A − 1, if 4 divides M then 4 divides A − 1
As a consequence, if M is square free, A = 1 is the only value, which leads to maximal period. If M = 3×2^k with k > 1, then A = 1 + 12n, with some integer n.--LaszloHars 22:18, 5 July 2007 (UTC)
- The case B=0 belongs to a different article now: Park-Miller RNG. Maxal (talk) 22:58, 22 February 2008 (UTC)
[edit] RNG from Numerical Recipes
...it is not generally realised that the above generator popularised by Numerical Recipes produces alternately odd and even results!
I have generated a small sample of numbers with this algorithm in python:
x=123492981749283749834750984327 for j in range(20): x=(1664525*x+1013904223)%(2E32) print x/2E32
which has an output:
0.785777231133 0.845651109565 0.413148398268 0.837627646619 0.158487884297 0.045609648011 0.399356437255 0.773722523776 0.483888193664 0.995558579190 0.144025803070 0.549855381779 0.029355124381 0.338410701578 0.073043339593 0.464837504736 0.647570258534 0.884586895494 0.002221975800 0.534268735330
which doesn't seem to alternate odd an even results. Maybe it's a problem in my implementation, can someone clarify this?
- Your code is incorrect. It should read:
x=123492981749283749834750984327 for j in range(20): x=(1664525*x+1013904223)%(2**32) print x
- 72.57.120.3 07:27, 16 July 2006 (UTC)
- 1. The initial value for x should be mod(123492981749283749834750984327,2^32) = 1746698375. It is unclear, how a larger than 64-bit number is handled.
2. 2E32 means 2*10^32, we need here 2^32 (2**32) 3. The result of the division is double precision floating point of 53 bits precision, which is converted to decimal for printing. The conversion hides the periods in the last bits, but they are still there. --LaszloHars 22:50, 5 July 2007 (UTC)
I reverted the removal of #'s 4 and 5 for exhibiting a full period. While they may be obvious to someone who is familiar to the problem, they may not be obvious to someone coming across LCG for the first time. Also, as far as I can tell, points 4 and 5 are not explicitly covered in points 1-3, and every time I've seen this list of 5 criteria in print, all 5 criteria are always listed. Halcyonhazard 18:18, 3 February 2007 (UTC)
[edit] LCRGs can be reasonable
I feel this article has a fairly strong anti-LCRG POV. For example, the problems with low order bits having short periods is only a problem with M is a power of two. Lets look at some C code which I just wrote (and dontate to the public domain):
/* Public domain */ /* These numbers are from page 6 of "Tables of Linear Congruential * Generators of Different Sizes and Good Lattice Structure" by Pierre * L'Ecuyer */ /* This is for a unsigned 32 bit number */ unsigned int gen(unsigned int a) { long long z; z = a; z *= 279470273; z %= 4294967291U; a = z; return a; } main() { unsigned int a; /* 32 bit unsigned */ int b; a = 1; for(b=0;b<100;b++) { a = gen(a); printf("%08x\n",a); } exit(0); }
Now, the numbers generated with this code do not have short periods in the low order bits. The main disadvantage of this code is that we have to use a temporary 64-bit number (the "long long" above) in order to generate a 32-bit random number, since we are doing a multiplication modulo an almost 32-bit number. I will modify the article accordingly; feel free to edit my edits. :) Samboy 17:28, 23 May 2007 (UTC)
[edit] better and faster pseudorandom number generators
This is a maximal period generator with the prime modulus 4294967291 = 2^32-5. It works, but it is still a very poor generator. Although the modulo operation can be performed without multiplication or division (with a shift-add, and compare), the overall work is more than with an Ax+B mod 2^32 type generator, and six 32-bit numbers are never generated (while the Ax+B generator has no missing values). The regularity is not better, only different.
If you want to break the simple patterns of the low order bits, you can drop a few entries, dependent on some bits, which are discarded by the modulo operator. I have been using this technique for more than 30 years (see, e.g. in http://www.autohotkey.com/forum/viewtopic.php?p=132576#132576). The result is comparable (in speed and randomness) to Marsaglia's MWC generator, which uses the previous carry to whiten the low order bits, but uses no static memory.
There are better and faster pseudorandom number generators, mainly designed for embedded applications. See: L. Hars and G. Petruska: Pseudorandom Recursions - Small and Fast Pseudorandom Number Generators for Embedded Applications. EURASIP Journal on Embedded Systems, vol. 2007, Article ID 98417, 13 pages, 2007. doi:10.1155/2007/98417. http://www.hindawi.com/GetArticle.aspx?doi=10.1155/2007/98417 --LaszloHars 23:21, 5 July 2007 (UTC)
[edit] RANDU link
Someone should do something about the RANDU link -- Annon —Preceding unsigned comment added by 72.85.32.253 (talk) 04:47, 17 January 2008 (UTC)