Fortran - Knuth: Linear Congruential Method

Author

Knuth: Linear Congruential Method

allelopath

2006-12-11, 7:12 pm

I'm trying to figure out, possibly update, some old, old fortran code.
The function is a random number generator.
It has this one comment:
c i changed the values that knuth suggested to allow
c running on a 16 bit integer machine
Looking in The Art of Computer Programming Volume by Donald E. Knuth,
I guess this function uses the Linear Congruential Method described
there:
x(n+1) = (a * x(n) + c) mod m, n >= 0
where:
m = modulus, 0 < m
a = multplier, 0 <= a < m
c = increment, 0 <= c < m
x0 = starting value, 0 <= x0 < m

I can't figure out what values he changed.
I would assume the modulus at the very least.
I need values for 32 bit.

Can anyone offer advice?

------------------------------------------------------------------

integer itab(64)
logical init
integer am,ixold,iyold,ictr,ixoldi,iyoldi,ia,ic,
m
integer k,i
save itab, init, ictr
data init/.false./
data ixoldi/57721/
data iyoldi/17810/
data m/566927/
data ia/3612/
data ic/5701/

integer mod

c ... The table, and the value for init need to be saved.
if ((.not. init) .or. k.lt.0) then
am = m
init = .true.
ixold = ixoldi
iyold = iyoldi
do 99 i=1,64
ixold = mod((ia*ixold + ic),m)
itab(i) = ixold
99 continue
endif
if (k.le.0) ictr = 0
ictr = ictr + 1
k = ictr
iyold = mod((ia*iyold + ic),m)
ixold = mod((ia*ixold + ic),m)
i = 1 + mod(iyold,64)
randu2 = float(itab(i)) / am
itab(i) = ixold

return
end

Terence

2006-12-13, 4:18 pm

No programmed "random" number generator is really random.
And any set of such numbers cannot be guaranteed to cover the requested
range without some values never occurring.

I remember using amplified shot noise from a thermal resistor as a form
of random input to a physics experiment for estimating physical
constants (Boltzmann).
Maybe that's why you can legally get hold of small quantities of
radioactive isotopes; another source of random interval events that
follow a Gausian bell curve histogram..

I coded Knuth's method in assembler to call it as a library function,
because there was a lot of bit section extraction from the middle of
long integers after a series of mutiplications. It was easier and
faster then Fortran at the time.

e p chandler

2006-12-13, 4:18 pm

allelopath wrote:
> I'm trying to figure out, possibly update, some old, old fortran code.
> The function is a random number generator.
> It has this one comment:
> c i changed the values that knuth suggested to allow
> c running on a 16 bit integer machine

Ignore the comment. The code below runs in 31 bits.

> Looking in The Art of Computer Programming Volume by Donald E. Knuth,
> I guess this function uses the Linear Congruential Method described
> there:
> x(n+1) = (a * x(n) + c) mod m, n >= 0
> where:
> m = modulus, 0 < m
> a = multplier, 0 <= a < m
> c = increment, 0 <= c < m
> x0 = starting value, 0 <= x0 < m
>
> I can't figure out what values he changed.
> I would assume the modulus at the very least.
> I need values for 32 bit.

You mean Knuth, vol 2, page 32?

> Can anyone offer advice?

Yes, avoid using random number routines of unknown origin for anything
more sophisticated than moving the fish around in a screen saver.

>
> ------------------------------------------------------------------
>
> integer itab(64)
> logical init
> integer am,ixold,iyold,ictr,ixoldi,iyoldi,ia,ic,
m
> integer k,i
> save itab, init, ictr
> data init/.false./
> data ixoldi/57721/
> data iyoldi/17810/
> data m/566927/
> data ia/3612/
> data ic/5701/
>
> integer mod
>
> c ... The table, and the value for init need to be saved.
> if ((.not. init) .or. k.lt.0) then
> am = m
> init = .true.
> ixold = ixoldi
> iyold = iyoldi
> do 99 i=1,64
> ixold = mod((ia*ixold + ic),m)
> itab(i) = ixold
> 99 continue
> endif
> if (k.le.0) ictr = 0
> ictr = ictr + 1
> k = ictr
> iyold = mod((ia*iyold + ic),m)
> ixold = mod((ia*ixold + ic),m)
> i = 1 + mod(iyold,64)
> randu2 = float(itab(i)) / am
> itab(i) = ixold
>
> return
> end

Does any one have a source or reference for this code?

-- elliot

allelopath

2006-12-13, 4:18 pm

>>You mean Knuth, vol 2, page 32?
vol 2, *section* 3.2.1

e p chandler

2006-12-13, 4:18 pm

On Dec 11, 10:24 pm, "allelopath" <jda...@gmail.com> wrote:[color=darkred]

Section 3.2.1 does describe the general principles of the linear
congruential method but section 3.2.2 (page 32) refers to algorithm M,
and algorithm B, "randomizing by shuffling".

I have seen this code posted in one example in some code from netlib.
It also appeared in a few game programs during a web search. I'm not
sure that the code is known to be bad, after all the other constants
*are* relatively prime to the modulus.... but there are numerous other
well tested and well documented routines abounding in the literature
and on the net. If small word size is a constraint there are combined
generators form L'Ecuyer or Wichmann and Hill.

allelopath

2006-12-13, 4:18 pm

These numbers:
data ixoldi/57721/
data iyoldi/17810/
data m/566927/
data ia/3612/
data ic/5701/

57721 = 197.293
17810 = 2.5.13.137
566927 = 23.157.157
3612 = 2.2.3.7.43
5701 is the only prime

Where do they come from?
How do you know they are for use with 32 bit?

e p chandler

2006-12-13, 4:18 pm

[regarding a lcg random number program]

On Dec 12, 3:31 am, "allelopath" <jda...@gmail.com> wrote:
> These numbers:
> data ixoldi/57721/
> data iyoldi/17810/
> data m/566927/
> data ia/3612/
> data ic/5701/
>
> 57721 = 197.293
> 17810 = 2.5.13.137
> 566927 = 23.157.157
> 3612 = 2.2.3.7.43
> 5701 is the only prime

The good news is that none of them have any factors in common with the
modulus m. The bad news is that the modulus is small. As a result, the
period of the sequence produced (without shuffling) is relatively short
and sets (say triplets) of points are not too well dispersed in space.

>
> Where do they come from?

It's magic!

> How do you know they are for use with 32 bit?

Residues mod m run from 0 to m-1. The maximum value produced inside the
LCG part of the routine is (m-1)*ia + ic which is less than 2**31 (2 to
the 31st power) So this code would run in 32 bits.

Gordon Sande

2006-12-13, 4:18 pm

On 2006-12-11 21:52:15 -0400, "Terence" <tbwright@cantv.net> said:

> No programmed "random" number generator is really random.

That is why the usual/correct name is PSEUDO random number generator.

> And any set of such numbers cannot be guaranteed to cover the requested
> range without some values never occurring.

The usual PRNGs can be analyzed with suitable mathematics. Typically
quite clasical number theory which many folks find dificult. One of
the simpler results is whether the PRNG in fact generates all possible
values. In one dimension the result is trivial to the extent that no
even vaguely plausible PRNG violates this. The property is called
one dimensional equidistribution.

In two dimensions the question is whether all possible pairs are generated.
If the PRNG is a LINEAR congruence the answer is clearly no. But if it
is a second order congruence the result will be true for even vaguely
plausible PRNGs. The question of how good a linear recurence can be for
two dimensions is addressed by the number theory topic called "geometry
of numbers" and leads to the full cycle test like the "Specral Test".
If the values in the pairs are of lower precision then the result can
be true. The number theory tells us how much precision is possible before
pairs start to be missing. Poor PRNGs can fall rather far short of the
possiblities. And so on for higher orders and dimensions ...

Many seemingly simple simulations have a requirement for precision and
dimension that can not be meet with typical simple linear congruence
based PRNGs. There are a number of very long cycle PRNGs that can be
used.

It is not clear what you were trying to say as the result is either:
- a trivial result of counting in simple cases
- a misleading comment, even false in many cases, on the issue of quality
of PRNGs.

> I remember using amplified shot noise from a thermal resistor as a form
> of random input to a physics experiment for estimating physical
> constants (Boltzmann).

Instrumentation for physical based random numbers is often a subtle and
dificult problem. Lots of documented problems.

> Maybe that's why you can legally get hold of small quantities of
> radioactive isotopes; another source of random interval events that
> follow a Gausian bell curve histogram..

How do you tell if the event is a positive interval after the previous
one or a negative interval before it? Remember the Gaussian has infinite
range from very negative to very posive.

> I coded Knuth's method in assembler to call it as a library function,
> because there was a lot of bit section extraction from the middle of
> long integers after a series of mutiplications. It was easier and
> faster then Fortran at the time.

Knuth is an excellent introduction that is now 25 years old and a bit dated
for some classes of results. Knuth offer several families of methods so
it is unclear what the "Knuth method" might be.

allelopath

2006-12-13, 4:18 pm

getting further into it:

The variable "am" seems to be superfluous.
It looks as though it m can be used instead.
or am I missing something?

the array size, referred to as k in the book, is "some number chosen
for convenience, usually in the neighborhood of 100"
This code uses 64, which depending on your point of view may or may not
be convenient and may or may not be in the neighborhood of 100.
64 is misleading, i think, as one might try to link it to 64 bit
processing or 2**6, but I don't think it is related to these.
or again, am I missing something?

glen herrmannsfeldt

2006-12-13, 4:18 pm

allelopath <jdaues@gmail.com> wrote:
> getting further into it:
(previously snipped PRNG algorithm)

> the array size, referred to as k in the book, is "some number chosen
> for convenience, usually in the neighborhood of 100"

This is the equivalent of using multiple decks of cards to
reduce the non-randomness of each individual deck.

> This code uses 64, which depending on your point of view may or may not
> be convenient and may or may not be in the neighborhood of 100.
> 64 is misleading, i think, as one might try to link it to 64 bit
> processing or 2**6, but I don't think it is related to these.
> or again, am I missing something?

Yes, the 64 is unrelated to 64 bit processors. There are
64 PRNGs, one is chose at random (using another PRNG), and then
the next value of the chosen one is computed.

-- glen

e p chandler

2006-12-13, 4:18 pm

On Dec 12, 2:50 pm, "allelopath" <jda...@gmail.com> wrote:
> getting further into it:
>
> The variable "am" seems to be superfluous.
> It looks as though it m can be used instead.
> or am I missing something?
>
> the array size, referred to as k in the book, is "some number chosen
> for convenience, usually in the neighborhood of 100"
> This code uses 64, which depending on your point of view may or may not
> be convenient and may or may not be in the neighborhood of 100.
> 64 is misleading, i think, as one might try to link it to 64 bit
> processing or 2**6, but I don't think it is related to these.
> or again, am I missing something?

64 is chosen as a convenient power of two. This makes it easy to
extract the upper bits of a binary "random" number and stuff them into
an index register. The remaining lower bits can then be used to replace
the value from the bucket you are currently emptying as addressed by
that register. (Bearing in mind that the lowest bits of *some* PRNGs
cycle regularly.)

allelopath

2006-12-13, 4:18 pm

Thanks for all the help, e p chandler and to all the others too. Smart
group of folks here.

allelopath

2006-12-15, 7:06 pm

I've rewritten the original function.
Removed some extraneous statements.
Renamed the variables to something more coherent.
The original code was probably f77, this new code is not.
Any critique is welcome. Thanks.

----------------------------------------------------------------
real function randu2(k)

c ****************************************
************

c modulus function
integer mod

integer k
integer i
integer am

integer, parameter :: START_X = 57721
integer, parameter :: START_Y = 17810
integer, parameter :: MODULUS = 566927
integer, parameter :: MULTIPLIER = 3612
integer, parameter :: INCREMENT = 5701
integer, parameter :: ARRAY_SIZE = 64

integer, save :: v_table(ARRAY_SIZE)
integer, save :: counter
integer, save :: x_sequence
integer, save :: y_sequence

logical, save :: initialized=.false.

c ... The table, and the value for 'initialized' need to be saved.
c enter this if body once per run, on first call
if ((.not. initialized) .or. k.lt.0) then
initialized = .true.
x_sequence = START_X
y_sequence = START_Y

do 99 i=1,ARRAY_SIZE
x_sequence =
- mod((MULTIPLIER * x_sequence + INCREMENT),MODULUS)
v_table(i) = x_sequence
99 continue
endif

counter = counter + 1
k = counter

c generate next numbers in sequence
y_sequence = mod((MULTIPLIER * y_sequence + INCREMENT),MODULUS)
x_sequence = mod((MULTIPLIER * x_sequence + INCREMENT),MODULUS)

i = 1 + mod(y_sequence,ARRAY_SIZE)

c generate the random number (return this value)
randu2 = float(v_table(i)) / MODULUS

v_table(i) = x_sequence

c write (*,*) randu2

return
end

allelopath

2006-12-15, 7:06 pm

I should also statement that it does produce seemingly random numbers.
Also, i see that I forgot to remove the declaration of the integer am

allelopath

2006-12-15, 7:06 pm

Craig Powers

2006-12-15, 7:06 pm

allelopath wrote:
> I should also statement that it does produce seemingly random numbers.
> Also, i see that I forgot to remove the declaration of the integer am

How have you evaluated "seemingly"? e.g. the naked eye is a poor way to
check for issues of bad randomness.

allelopath

2006-12-15, 7:06 pm

I plotted 5000 values on a graph.
I'm not saying that this is scientifically rigorous.

On Dec 15, 12:41 pm, Craig Powers <eni...@hal-pc.org> wrote:
> allelopath wrote:
> check for issues of bad randomness.

robin

2006-12-15, 7:06 pm

allelopath wrote in message <1165871974.073594.37380@j44g2000cwa.googlegroups.com>...
>I'm trying to figure out, possibly update, some old, old fortran code.
>The function is a random number generator.
>It has this one comment:
>c i changed the values that knuth suggested to allow
>c running on a 16 bit integer machine
>Looking in The Art of Computer Programming Volume by Donald E. Knuth,
>I guess this function uses the Linear Congruential Method described
>there:
>x(n+1) = (a * x(n) + c) mod m, n >= 0
> where:
> m = modulus, 0 < m
> a = multplier, 0 <= a < m
> c = increment, 0 <= c < m
> x0 = starting value, 0 <= x0 < m
>
>I can't figure out what values he changed.

Why not compare the constants in the code in Knuth's book with those
in the 16-bit code?

>I would assume the modulus at the very least.
>I need values for 32 bit.
>
>Can anyone offer advice?
>
>------------------------------------------------------------------
>
> integer itab(64)
> logical init
> integer am,ixold,iyold,ictr,ixoldi,iyoldi,ia,ic,
m
> integer k,i
> save itab, init, ictr
> data init/.false./
> data ixoldi/57721/
> data iyoldi/17810/
> data m/566927/
> data ia/3612/
> data ic/5701/

allelopath

2006-12-15, 7:06 pm

>> the 16-bit code?
if you dig through this thread, e p chandler showed that the code was
in reality written for 32-bit, despite the 16-bit comment.

In my code critique request, I guess I was looking for comments
regarding fortran conventions or perhaps reasonable optimizations.

Craig Powers

2006-12-18, 7:08 pm

allelopath wrote:
> I plotted 5000 values on a graph.
> I'm not saying that this is scientifically rigorous.
>
> On Dec 15, 12:41 pm, Craig Powers <eni...@hal-pc.org> wrote:
>

Nor even particularly non-scientifically rigorous. :)

I suppose it would catch really blindingly obvious problems, but there
are relatively common issues (e.g. poor randomness in the least
significant bit) that I'd guess it would miss.

Sponsored Links