C PCG64

In this appendix we describe our focal random number generator, PCG64, in more detail. A background discussion of binary integer representation and the type of bitwise operators used in the generator is provided in the previous appendix.

The permuted congruential random number generator (PCG) concept was pioneered by O’Neill (2014). The source code for various members of the family is hosted on GitHub at https://github.com/imneme/pcg-python.

PCG64 is a particular PCG choice that maintains a 128-bit internal state while producing 64-bit outputs. Obtaining a value from the generator requires two steps. First, the state is advanced with a linear congruential generator having multiplier 47026247687942121848144207491837523525 and an odd increment. Then, the final output is produced from using an xorshift and a rotation.

The LCG update phase serves the purpose of giving the generator a long period and ensuring that all states are visited. The xorshift then provides a quick, non-linear, “avalanche” that mixes high- and low-order bits. The subsequent rotation prevents fixed correlations between bit positions and spreads entropy so that low-order output bits are just as random as the high-order bits.

When taken together, all the transformations in PCG64 yield a generator with a simple state transition and a high-quality 64-bit output that is known to pass all the TestU01 performance criteria.

C.1 A Toy Generator

To illustrate the basic structure of PCGs we will explore a small example having an 8-bit state and a 4-bit output. While this generator is not suitable for real simulation work, it allows us to see the two key components of PCGs at work.

The core of any PCG is its congruential component. For our 8-bit example we will use \[\begin{equation} state_{t+1} = (45 \ast state_t + 77 ) \bmod 2^8 \label{eq:pcg8-lcg} \end{equation}\] The generator constants \(a\) = 45 and \(c\) = 77 are somewhat arbitrary; but, the Hull-Dobell theorem of the previous appendix ensures that the generator has the full 256 period: i.e., \(c\) = 77 is odd, so \(m = 2^8\) and \(c\) are coprime, and \(a - 1\) = 44 is divisible by 4 (and hence by 2, the only prime factor of \(m\)).

Python code that implements the generator looks like

def xorshifted(x, a, b):
    return ((x >> a) ^ x) >> b

def rotate_right(value, rot):
    # Rotate an 8-bit value right by at most 7 bits
    rot %= 8
    return ((value >> rot) | ((value << (8 - rot)) & 255))

def pcg8(state):

    # Linear congruential update
    new_state = (state * 45 + 77) % 256

    # Mix bits with an xorshift
    mixed = xorshifted(new_state, 2, 3)

    # Derive a rotation amount from the top 3 bits of the new state
    rot = new_state >> 5  # Value between 0 and 7

    # Rotate the mixed value masked to one byte
    output = rotate_right(mixed & 255, rot)

    # return the congruential update and generator output
    return new_state, output

Here is how the toy generator works. First the generator’s internal state is updated using the linear congruential formula. The new state is given by

\[\text{new\_state} = (\text{old\_state} \ast 45 + 77)\;\bmod 256.\]

The new state is XORed with itself, shifted right by two bits and the result of that is shifted right by three more bits. \[ \text{mixed} = \bigl((\text{new\_state} >> 2) \text{\textasciicircum} \text{new\_state}\bigr) >> 3 \]

This serves the purpose of mixing together bits from different positions.

The three most significant bits of the new state determine how many bits to rotate the mixed value to the right. A rotation moves low-order bits around to high positions and vice versa, further scrambling the output.

Let us work through a numerical example in Python to illustrate all the generator’s steps. For that purpose take the input state to be 42. We will compute two successive outputs to illustrate the internal transformations. Recall that all arithmetic is performed modulo 256 and that shifts and rotations operate on 8-bit values.

For the first iteration the linear congruential update is

new_state = (42*45 + 77) % 256
new_state

## 175

This is subjected to XOR–shift mixing to obtain

mixed = ((175 >> 2) ^ 175) >> 3
mixed

## 16

Now,

bin(175)

## '0b10101111'

So, the top three bits of 175 are 101 (5 in decimal) which gives 5 as the rotation amount and

output = (16 >> 5) | (16 << 3)
output

## 128

is the random number returned by the generator.

For the second iteration we feed the new state 175 into the generator to get

new_state = (175*45 + 77) % 256
new_state

## 16

as the congruential update. Now xorshift to obtain

mixed = ((16 >> 2) ^ 16) >> 3
mixed

## 2

Since, 16 is 00010000 in binary, the rotation amount is zero which leaves mixed unchanged. Consequently, the generator returns 2 as its next random number.

Comparing the LCG generator states with the generator’s output gives

bin(175)
bin(128)
bin(16)
bin(2)

## '0b10101111'
## '0b10000000'
## '0b10000'
## '0b10'

These two iterations of the toy generator demonstrate how the state changes and how the xorshift and rotation steps mix the bits to produce outputs that differ markedly from the state itself. Continuing in this fashion produces a sequence that cycles through all 256 states before repeating.

Because the state space has only 256 values, we can enumerate the entire cycle of our simple generator. Figure C.1 shows the output from PCG8 as well as the original congruential generator. The permuted output shows less structure than for the LCG.

Figure C.1: Congruential vs permuted generator

It is not difficult now to understand the workings of the full PCG64 generator. Just like PCG8, this generator advances the state using a congruential generator. Then, it XOR shifts it and rotates the result by an amount determined by the previous LCG state.

C.2 Independent Streams

There are two reasons we have focused on PCG64 in this book: 1) it has an excellent performance record as a uniform number generator and 2) it is easy to use in parallel processing.

For a generator to be effective in parallel we need the generated number streams to be distinct across processors. PCG generators accomplish this. The increments of their LCGs act as stream selectors. Each odd increment gives a full-period sequence and, as noted in the previous appendix, changing the increment alters the underlying orbit of states thereby producing distinct streams. These streams are deterministic and therefore not “independent” in any probabilistic sense. But, they are constructed to be distinct sequences and to have good empirical statistical behavior.

Each admissible odd 128-bit increment selects a different transition rule. The resulting full-period state sequences may visit the same individual states, but sequences with different increments cannot be shifted copies of one another. PCG64 has period \(2^{128}\) for each stream and supports \(2^{127}\) possible streams. These streams are deterministic; “independence” refers to their construction and empirical statistical behavior, not probabilistic independence. Statistical tests of this between-stream behavior are developed in Ismay (2013).

C.3 PCG64 in R and Python

In R an implementation of PCG64 is provided in the dqrng package. To use it we need to first load the package

library(dqrng)

and specify the random number generator kind via

dqRNGkind("pcg64")

The seed for PCG64 is a vector of two 32-bit integers that can be obtained using the generateSeedVectors function from the dqrng package. In general, generateSeedVectors(n_seed) will produce n_seed integer vector pairs. Of course, generateSeedVectors requires a random seed in its own right and that is obtained from R’s native random number generator. So, for example,

set.seed(1234)
generateSeedVectors(1)

## [[1]]
## [1] -1622172679  -597384213

generateSeedVectors(1)

## [[1]]
## [1]   998812292 -1315992701

set.seed(1234)
generateSeedVectors(1)

## [[1]]
## [1] -1622172679  -597384213

There are a couple of things this output demonstrates. First, generateSeedVectors returns the 2-dimensional vector of integers for seeding the generator in a list. Secondly, and somewhat surprisingly, it can produce negative numbers. But, the seeds for the generator are supposed to be non-negative!

The “problem” here is all related to how one interprets the bits for an integer’s representation. R stores integers in a 32-bit 2’s complement format with the most significant bit determining the sign. In contrast, the C++ implementation that underlies the code in dqrng views values received from R as unsigned 32-bit integers. The sign bit simply becomes the coefficient for \(2^{31}\) in their binary representation.

Setting the seed and increment for PCG64 in R requires two 2-dimensional vectors of 32-bit integers: one vector for the seed and one for the increment. We have already seen that generateSeedVectors can be used to provide such quantities. The vectors obtained in this manner will suffice for setting the seed of the generator and they also can be used to set the increment or stream. There is a safeguard in the underlying C++ dqrng code that ensures the increment that is used is odd.

PCG64 is the default random number generator in Python. The seed and increment are derived using SeedSequence, which uses a hashing technique to convert any user input into a high-quality 128-bit internal state and a fixed odd 128-bit increment. For parallel work the PCG64 bit generator utilizes the NumPy.random SeedSequence.spawn method. This is the strategy employed in sim_template.py in Chapter 4.

B Random Number Generators

References