The Random Number Generator

In recent years, progress has been made in constructing portable generators with large periods and other good properties; see the review [Jam90]. Therefore the current version contains a random number generator based on the algorithm proposed by Marsaglia, Zaman and Tsang [Mar90]. This routine should work on any machine with a mantissa of at least 48 digits, i.e. on computers with a 64-bit (or more) representation of double precision real numbers. Given the same initial state, the sequence will also be identical on different platforms. This need not mean that the same sequence of events will be generated, since the different treatments of roundoff errors in numerical operations will lead to slightly different real numbers being tested against these random numbers in IF statements. Also code optimization may lead to a divergence of the event sequence.

Apart from nomenclature issues, the coding of PYR as a function rather than a subroutine, and the extension to double precision, the only difference between our code and the code given in [Jam90] is that slightly different algorithms are used to ensure that the random number is not equal to 0 or 1 within the machine precision. Further developments of the algorithm has been proposed [Lüs94] to remove residual possibilities of small long-range correlations, at the price of a slower generation procedure. However, given that PYTHIA is using random numbers for so many different tasks, without any fixed cycle, this has been deemed unnecessary.

The generator has a period of over $10^{43}$, and the possibility to obtain almost $10^9$ different and disjoint subsequences, selected by giving an initial integer number. The price to be paid for the long period is that the state of the generator at a given moment cannot be described by a single integer, but requires about 100 words. Some of these are real numbers, and are thus not correctly represented in decimal form. The old-style procedure, which made it possible to restart the generation from a seed value written to the run output, is therefore not convenient. The CERN library implementation keeps track of the number of random numbers generated since the start. With this value saved, in a subsequent run the random generator can be asked to skip ahead the corresponding number of random numbers. PYTHIA is a heavy user of random numbers, however: typically 30% of the full run time is spent on random number generation. Of this, half is overhead coming from the function call administration, but the other half is truly related to the speed of the algorithm. Therefore a skipping ahead would take place with 15% of the time cost of the original run, i.e. an uncomfortably high figure.

Instead a different solution is chosen here. Two special routines are provided for writing and reading the state of the random number generator (plus some initialization information) on a sequential file, in a platform-dependent internal representation. The file used for this purpose has to be specified by you, and opened for read and write. A state is written as a single record, in free format. It is possible to write an arbitrary number of states on a file, and a record can be overwritten, if so desired. The event generation loop might then look something like:

28.

save the state of the generator on file (using flag set in point 3 below),

29.

generate an event,

30.

study the event for errors or other reasons why to regenerate it later; set flag to overwrite previous generator state if no errors, otherwise set flag to create new record;

31.

loop back to point 1.

With this procedure, the file will contain the state before each of the problematical events. These events can therefore be generated in a shorter run, where further information can be printed. (Inside PYTHIA, some initialization may take place in connection with the very first event generated in a run, so it may be necessary to generate one ordinary event before reading in a saved state to generate the interesting events.) An alternative approach might be to save the state every 100 events or so. If the events are subsequently processed through a detector simulation, you may have to save also other sets of seeds, naturally.

Unfortunately, the procedure is not always going to work. For instance, if cross section maximum violations have occured before the interesting event in the original run, there is a possibility that another event is picked in the re-started one, where the maximum weight estimate has not been updated. Another problem is the multiple interaction machinery, where some of the options contain an element of learning, which again means that the event sequence may be broken.

In addition to the service routines, the common block which contains the state of the generator is available for manipulation, if you so desire. In particular, the initial seed value is by default 19780503, i.e. different from the Marsaglia/CERN default 54217137. It is possible to change this value before any random numbers have been generated, or to force re-initialization in mid-run with any desired new seed.

It should be noted that, of course, the appearance of a random number generator package inside PYTHIA does in no way preclude the use of other routines. You can easily revert to having PYR as nothing but an interface to an arbitrary external random number generator; e.g. to call a routine RNDM all you need to have is

      FUNCTION PYR(IDUMMY)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
  100 PYR=RNDM(IDUMMY)
      IF(PYR.LE.0D0.OR.PYR.GE.1D0) GOTO 100
      RETURN
      END

The random generator subpackage consists of the following components.

Purpose:

to generate a (pseudo)random number R uniformly in the range 0$<$R$<$1, i.e. excluding the endpoints.

IDUMMY :

dummy input argument; normally 0.

Purpose:

to dump the current state of the random number generator on a separate file, using internal representation for real and integer numbers. To be precise, the full contents of the PYDATR common block are written on the file, with the exception of MRPY(6).

LFN :

(logical file number) the file number to which the state is dumped. You must associate this number with a true file (with a platform-dependent name), and see to it that this file is open for write.

MOVE :

choice of adding a new record to the file or overwriting old record(s). Normally only options 0 or $-$1 should be used.

= 0 (or > 0) :

add a new record to the end of the file.

= -1 :

overwrite the last record with a new one (i.e. do one BACKSPACE before the new write).

= $-n$ :

back up $n$ records before writing the new record. The records following after the new one are lost, i.e. the last $n$ old records are lost and one new added.

Purpose:

to read in a state for the random number generator, from which the subsequent generation can proceed. The state must previously have been saved by a PYRGET call. Again the full contents of the PYDATR common block are read, with the exception of MRPY(6).

LFN :

(logical file number) the file number from which the state is read. You must associate this number with a true file previously written with a PYRGET call, and see to it that this file is open for read.

MOVE :

positioning in file before a record is read. With zero value, records are read one after the other for each new call, while non-zero values may be used to navigate back and forth, and e.g. return to the same initial state several times.

= 0 :

read the next record.

= $+n$ :

skip ahead $n$ records before reading the record that sets the state of the random number generator.

= $-n$ :

back up $n$ records before reading the record that sets the state of the random number generator.

Purpose:

to contain the state of the random number generator at any moment (for communication between PYR, PYRGET and PYRSET), and also to provide you with the possibility to initialize different random number sequences, and to know how many numbers have been generated.

MRPY(1) :

(D=19780503) the integer number that specifies which of the possible subsequences will be initialized in the next PYR call for which MRPY(2)=0. Allowed values are 0$\leq$MRPY(1)$\leq$900000000, the original Marsaglia (and CERN library) seed is 54217137. The MRPY(1) value is not changed by any of the PYTHIA routines.

MRPY(2) :

(D=0) initialization flag, put to 1 in the first PYR call of run. A re-initialization of the random number generator can be made in mid-run by resetting MRPY(2) to 0 by hand. In addition, any time the counter MRPY(3) reaches 1000000000, it is reset to 0 and MRPY(2) is increased by 1.

MRPY(3) :

(R) counter for the number of random numbers generated from the beginning of the run. To avoid overflow when very many numbers are generated, MRPY(2) is used as described above.

MRPY(4), MRPY(5) :

I97 and J97 of the CERN library implementation; part of the state of the generator.

MRPY(6) :

(R) current position, i.e. how many records after beginning, in the file; used by PYRGET and PYRSET.

RRPY(1) - RRPY(97) :

the U array of the CERN library implementation; part of the state of the generator.

RRPY(98) - RRPY(100) :

C, CD and CM of the CERN library implementation; the first part of the state of the generator, the latter two constants calculated at initialization.