Talk:Linear congruential generator

From Wikipedia, the free encyclopedia
Jump to: navigation, search
WikiProject Computer science (Rated Start-class, Mid-importance)
WikiProject icon This article is within the scope of WikiProject Computer science, a collaborative effort to improve the coverage of Computer science related articles on Wikipedia. If you would like to participate, please visit the project page, where you can join the discussion and see a list of open tasks.
Start-Class article Start  This article has been rated as Start-Class on the project's quality scale.
 Mid  This article has been rated as Mid-importance on the project's importance scale.
 

C code example of LCG for 32-bit CPU[edit]

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?

LCG examples[edit]

  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 2^{31}-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 2^{31}-250, which is also the maximum one for Park-Miller RNG modulo prime 2^{31} - 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)

Grammar issues[edit]

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)

disproof[edit]

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.


I tried using a=184,c=4 and m=549. In spite of the fact that all conditions are indeed satisfied, the sequence we get is not random when plotted against time or loop variable! In fact its too regular!! — Preceding unsigned comment added by 210.212.187.69 (talk) 12:40, 11 October 2011 (UTC)

Rules for the maximum period[edit]

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.

Case B=0[edit]

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)

RNG from Numerical Recipes[edit]

...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)

LCRGs can be reasonable[edit]

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)

better and faster pseudorandom number generators[edit]

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)


RANDU link[edit]

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)

output bits of seed in rand() / Random(L)[edit]

"output bits of seed in rand() / Random(L)" What is rand()? What is Random(X)? What is L? Why are these divided? What is meant by "output bits"? How does a seed have "output bits" if it's just a number? 130.89.1.11 (talk) 13:56, 25 June 2008 (UTC)

  • rand() is a standard C/C++ function, Random(L) is a standard Pascal/Delphi function (and L is its argument). Their implementations often output/return only part (e.g., low or high half) of the seed. Maxal (talk) 04:14, 9 February 2012 (UTC)

References and footnotes[edit]

I notice that the references and footnotes are together under the section heading "References". It is does not look well (from an aesthetic point of view). I will introduce a Notes section proper.--Михал Орела (talk) 06:43, 14 January 2009 (UTC)

good values for other modulus[edit]

Are there "proven" good values for a and c for a modulus m = 264, 2128 or 2256 and other "big" moduli? --82.83.126.246 (talk) 09:51, 18 May 2009 (UTC)

  • See conditions in "Period length" section. Maxal (talk) 04:15, 9 February 2012 (UTC)

ZX81/ZX Spectrum RND function[edit]

From the manual (the text is approximately the same in each version, with only the model of computer changed):

Let p be a (large) prime, and let a be a primitive root \pmod{p}.

Then if b_i is the residue of a^i \pmod{p} (1 \le b_i \le p-1), the sequence

\frac{b_i-1}{p-1}

is a cyclical sequence of p-1 distinct numbers in the range 0 to 1 (excluding 1). By choosing a suitably, these can be made to look fairly random.

65537 is a Fermat prime, 2^{16}+1. Because the multiplicative group of non-zero residues modulo 65537 has a power of 2 as its order, a residue is a primitive root if and only if it is not a quadratic residue. Use Gauss' law of quadratic reciprocity to show that 75 is a primitive root modulo 65537.

The ZX81/ZX Spectrum uses p=65537 and a=75, and stores some b_i-1 in memory. RND entails replacing b_i-1 in memory by b_{i+1}-1, and yielding the result \frac{b_{i+1}-1}{p-1}.

(Note that both computers actually perform the calculation in full 32-bit floating-point arithmetic rather than integer modulo arithmetic, which would have been reasonably straightforward for a modulus of 65537 even on an 8-bit processor such as the Z80 used in the ZX81/ZX Spectrum.)

80.175.250.218 (talk) 19:25, 7 August 2009 (UTC)

Equation symbols[edit]

Surely it should be X_{n+1} = \left( a X_n + c \right)... rather than X_{n+1} \equiv \left( a X_n + c \right)...  ?Iapetus (talk) 10:31, 11 June 2013 (UTC)