Path: cactus.org!milano!cs.utexas.edu!qt.cs.utexas.edu!yale.edu!newsserver.
+     jvnc.net!howland.reston.ans.net!paladin.american.edu!news.univie.ac.at!
+     hp4at!mcsun!uknet!comlab.ox.ac.uk!imc
From: imc@comlab.ox.ac.uk (Ian Collier)
Newsgroups: comp.programming

Subject: Re: Random number generator algorithms
Keywords: LCGs
Message-ID: <3167.imc@uk.ac.ox.prg>
Date: 25 Feb 93 13:43:13 GMT
References: <1993Feb19.195525.23105@linus.mitre.org>
Organization: Oxford University Computing Laboratory
Lines: 47
X-Local-Date: Thursday, 25th February 1993 at 1:42pm GMT
Originator: imc@msc6.comlab

It looks like our news server is a few days behind at the moment, but...

In article <1993Feb19.195525.23105@linus.mitre.org>, pate@severian.mitre.org (La ura Pate) wrote:
>                                    I've implemented a simple linear
>congruential generator, Z(i) = (a*Z(i-1) + c ) ( mod m ). Law&Kelton
>give several guidelines for choosing a, c, and m to ensure full
>period, and IID random numbers. I have two questions:

>1. This algorithm is fairly old, are there new and better algorithms
>that I should be looking at?

It's a fairly good method whether it's old or not; however if m contains a
high power of 2 then the least significant few bits can be rather less than
perfectly random.  Knuth does give a couple of other interesting
algorithms.  One of these which you might look at is
X(i) = {X(i-p) + X(i-q)} mod m, for particular values of p and q (they are
given, but I can't remember them).  Obviously you will need to seed the
generator initially with max(p,q) random values.

>2. If this is a reasonable algorithm, does anyone have any suggestions
>for choosing a and c when m is constrained to be 2^16 (I have to
>multiple another value by m and the result must be < INT_MAX). The
>examples in Law&Kelton and those referenced in Knuth all have much
>larger m's of 2^32 or 2^35. 

To disobey your request slightly... with m=65537 (this is prime, so the
sequence does not suffer the disadvantage mentioned above) it is common to
choose a=75 and c=0, giving a purely multiplicative generator.  This gives
two cycles, one of period 65536 and the other containing the single element
0.  Each member of the 65536-cycle will clearly fit into 16 bits if you
decrement it.  So the sequence looks like this:  X(0) is any 16-bit number,
and X(i) = [(75*(X(i-1)+1)) mod 65537] - 1.

Multiplying a 16-bit number by 65537 is easy, as you just store it into each
of the high and low halfwords, and it cannot overflow.  Reducing a 32-bit
number mod 65537 is also easy, as Knuth observes: you divide the product
into its two sets of 16 bits and subtract the higher from the lower.  If the
answer is negative, increment the 16-bit 2's complement result as if it were
a positive number.

As you know, any linear congruential generator gives a cycle of length at
most m (when m is prime the most is m-1 unless you take a=c=1).  If you
require a larger cycle, Knuth gives algorithms for "shuffling" the numbers
obtained from a generator (called algorithm M, I think).

Ian Collier
Ian.Collier@prg.ox.ac.uk | imc@ecs.ox.ac.uk