“Finally, we give you Knuth's suggestion for a portable routine,
which we have translated to the present conventions as RAN3. This
is not based on the linear congruential method at all, but rather
on a ** subtractive method**. One might hope that its weaknesses,
if any, are therefore of a highly different character from the
weaknesses, if any, of RAN1 above [not given]. If you ever suspect
trouble with one routine, it is a good idea to try the other in the
same application. RAN3 has one nice feature: if your machine is
poor on integer arithmetic, (i.e. is limited to 16-bit integers),
substitution of the two "commented" lines for the one directly
following them will render the routine entirely floating point.”

FUNCTION RAN3(IDUM)

Returns a uniform random deviate between 0.0 and 1.0. Set IDUM to any negative value

to initialize or reinitialize the sequence.

C IMPLICIT REAL*4(M)

C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5E-7)

PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.E-9)

According to Knuth, any large MBIG, and any smaller (but still large) MSEED can be

substituted for the above values.

DIMENSION MA(55) Save MA. This value is special and should not be modified; see Knuth

DATA IFF /0/

IF(IDUM.LT.0.OR.IFF.EQ.0)THEN Initialization

IFF=1

MJ=MSEED-IABS(IDUM) Initialize MA(55) using the seed IDUM and the large number MSEED.

MJ=MOD(MJ,MBIG)

MA(55)=MJ

MK=1

DO 11 I=1,54 Now initialize the rest of the table

II=MOD(21*I,55) in a slightly random order

MA(II)=MK with numbers that are not especially random

MK=MJ-MK

IF(MK.LT.MZ)MK=MK+MBIG

MJ=MA(II)

11 CONTINUE

DO 13 K=1,4 We randomize them by "warming up the generator"

DO 12 I=1,55

MA(I)=MA(I)-MA(1+MOD(I+30,55))

IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG

12 CONTINUE

13 CONTINUE

INEXT=0 Prepare indices for our first generated number

INEXTP=31 The constant 31 is special; see Knuth

IDUM=1

ENDIF

INEXT=INEXT+1 Here is where we start, except on initialization. Increment INEXT,

IF(INEXT.EQ.56)INEXT=1 wrapping around 56 to 1.

INEXTP=INEXTP+1 Ditto for INEXTP

IF(INEXTP.EQ.56)INEXTP=1

MJ=MA(INEXT)-MA(INEXTP) Now generate a new random number subtractively

IF(MJ.LT.MZ)MJ=MJ+MBIG Be sure that it is in range

MA(INEXT)=MJ and output the derived uniform deviate

RAN3=MJ*FAC

RETURN

END