diff --git a/docs/source/methods/random_numbers.rst b/docs/source/methods/random_numbers.rst index 0cefc9156bc..4376bbdb0c5 100644 --- a/docs/source/methods/random_numbers.rst +++ b/docs/source/methods/random_numbers.rst @@ -7,7 +7,7 @@ Random Number Generation In order to sample probability distributions, one must be able to produce random numbers. The standard technique to do this is to generate numbers on the interval :math:`[0,1)` from a deterministic sequence that has properties that -make it appear to be random, e.g. being uniformly distributed and not exhibiting +make it appear to be random, e.g., being uniformly distributed and not exhibiting correlation between successive terms. Since the numbers produced this way are not truly "random" in a strict sense, they are typically referred to as pseudorandom numbers, and the techniques used to generate them are pseudorandom @@ -15,6 +15,11 @@ number generators (PRNGs). Numbers sampled on the unit interval can then be transformed for the purpose of sampling other continuous or discrete probability distributions. +There are many different algorithms for pseudorandom number generation. OpenMC +currently uses `permuted congruential generator`_ (PCG), which builds on top of +the simpler linear congruential generator (LCG). Both algorithms are described +below. + ------------------------------ Linear Congruential Generators ------------------------------ @@ -37,8 +42,8 @@ be generated with a method chosen at random. Some theory should be used." Typically, :math:`M` is chosen to be a power of two as this enables :math:`x \mod M` to be performed using the bitwise AND operator with a bit mask. The constants for the linear congruential generator used by default in OpenMC are -:math:`g = 2806196910506780709`, :math:`c = 1`, and :math:`M = 2^{63}` (see -`L'Ecuyer`_). +:math:`g = 2806196910506780709`, :math:`c = 1`, and :math:`M = 2^{63}` (from +`L'Ecuyer `_). Skip-ahead Capability --------------------- @@ -50,7 +55,8 @@ want to skip ahead :math:`N` random numbers and :math:`N` is large, the cost of sampling :math:`N` random numbers to get to that position may be prohibitively expensive. Fortunately, algorithms have been developed that allow us to skip ahead in :math:`O(\log_2 N)` operations instead of :math:`O(N)`. One algorithm -to do so is described in a paper by Brown_. This algorithm relies on the following +to do so is described in a `paper by Brown +`_. This algorithm relies on the following relationship: .. math:: @@ -58,15 +64,26 @@ relationship: \xi_{i+k} = g^k \xi_i + c \frac{g^k - 1}{g - 1} \mod M -Note that equation :eq:`lcg-skipahead` has the same general form as equation :eq:`lcg`, so -the idea is to determine the new multiplicative and additive constants in -:math:`O(\log_2 N)` operations. +Note that equation :eq:`lcg-skipahead` has the same general form as equation +:eq:`lcg`, so the idea is to determine the new multiplicative and additive +constants in :math:`O(\log_2 N)` operations. -.. only:: html - .. rubric:: References +-------------------------------- +Permuted Congruential Generators +-------------------------------- +The `permuted congruential generator`_ (PCG) algorithm aims to improve upon the +LCG algorithm by permuting the output. The algorithm works on the basic +principle of first advancing the generator state using the LCG algorithm and +then applying a permutation function on the LCG state to obtain the output. This +results in increased statistical quality as measured by common statistical tests +while exhibiting a very small performance overhead relative to the LCG algorithm +and an equivalent memory footprint. For further details, see the original +technical report by `O'Neill +`_. OpenMC uses the +PCG-RXS-M-XS variant with a 64-bit state and 64-bit output. -.. _L'Ecuyer: https://doi.org/10.1090/S0025-5718-99-00996-5 -.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf .. _linear congruential generator: https://en.wikipedia.org/wiki/Linear_congruential_generator + +.. _permuted congruential generator: https://en.wikipedia.org/wiki/Permuted_congruential_generator diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index 581d696176f..ca4467719c4 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -17,7 +17,8 @@ constexpr uint64_t prn_stride {152917LL}; // stride between particles //============================================================================== // 64 bit implementation of the PCG-RXS-M-XS 64-bit state / 64-bit output -// geneator Adapted from: https://github.com/imneme/pcg-c +// geneator Adapted from: https://github.com/imneme/pcg-c, in particular +// https://github.com/imneme/pcg-c/blob/83252d9c23df9c82ecb42210afed61a7b42402d7/include/pcg_variants.h#L188-L192 // @techreport{oneill:pcg2014, // title = "PCG: A Family of Simple Fast Space-Efficient Statistically Good // Algorithms for Random Number Generation", author = "Melissa E. O'Neill",