Overview
Standard random number generators allow you get the current member of the sequence and update the position to the next. They do not allow you to walk backwards, skip around or simply ask: What’s the n^{th} member?
This public domain single header file library (prns.h) implements basic functions for random access reads of a pseudorandom sequence. It is specifically designed for this purpose and all functions are roughly the same complexity. Although it is relatively fast and has quite reasonable statistically quality, there are stateoftheart generators which are both faster and higher quality. So it is only of interest if you need a generator which can perform these uncommon operations in constant time. It additionally support a second independent sequence.
The structure is explained in: Implementation details
If seeking and skipping are of reduced interest we can change the base generation technique (SEE: An alternate solution) from that initially presented.
Documentation
The library is composed of one typedef prns_t
which holds the state data for a given generator and a handful of inline functions which provide basic operations.
Position in sequence manipulation
Returns the current position in the sequence.
Set the current position in the sequence to pos.
Moves the current position by the relative value offset.
Reading sequence members
Returns the n^{th} member of the sequence.
Returns the member of the sequence at the current position.
Returns the member of the sequence at the current position and increments.
Returns the member of the sequence at the current position and decrements.
Down sequence
In addition to the main random number sequence of the above functions the library supports a second independent sequence. Its state is initalized from the current position in the main sequence by:
the result of which is not a random number but a state value. To read from the sequence:
As a strawman example the following queries exactly one value from the main sequence and either increments or decrements the position depending on the boolean forward. It can then reads some unknown number of values from the down sequence to perform its function.
Utils
The bit finalizing function. (SEE: bit finalizing)
Implementation details and configuration options
The generator is structured like most modernish PRNG and hashing functions. One function updates the state and a second performs mixing before returning the final result.
State update
To make all position related operations constant time requires a simple closedform state update which is O(1). Using the members of a Weyl sequence ^{1} for state fills this requirement. If we call $ n $ the position in the sequence, $ s_n $ the state at that position, $ W $ the Weyl constant and $ k $ the signed offset we’re updating the position then state updates can be expressed as:
\[\begin{equation} \label{eq:statePos} s_{n+k}=s_{n}+k\thinspace W \end{equation}\]The natural value for $ s_0 $ is zero. Fixing $ n $ and/or $ k $ gives how state is related to position. Solving the closed form for the position $ n $:
\[\begin{equation} \label{eq:position} n=W^{1}\thinspace s_n \end{equation}\]where $ W^{1} $ is the modulo inverse of $ W $. So all state operations become modulo integer addition and/or multiplication.
The down sequence follows the same structure simply with different Weyl constant.
The constants for the main sequence are macro defines: PRNS_WEYL
defines the forward $ W $ and PRNS_WEYL_I
defines $ W^{1} $. The down sequence only support query and increment and its constant is defined by PRNS_WEYL_D
.
Bit finalizing
The state update is so simple that the burden of producing a sequence with statistical random properties is on the shoulders of the bit finalizing function. David Stafford posted a blog entry^{2} in which he presents fourteen variants of cascading two right xorshift multiples followed by an right xorshift:
and two of the variant have good performance for low entropy input mix01 and mix13 and a Weyl sequence is low entropy. A right xorshift can be expressed as matrix equation over $ \mathbb{F_2} $:
\[x' = (I+R^a)x\]where $ I $ is the identity and $ R $ the (right) shift matrices and $ a $ is the shift amount.
So the overall structure of a produced random number is alternating operations in $ \mathbb{Z_{2^{64}}} $ (Weyl sequence and products in the mixing function) and $ \mathbb{F_2} $.
The logical of this kind of construction is informally as follows. A modulo integer product with an odd constant produces no mixing on the lowest bit and the most on the highest. Each bit position from the lowesttohighest gets increasing more mixed. If this isn’t obvious then consider the long hand product in binary. Right xorshifts are typcially formed so the shift amount is around half the working bitwidth (~16 for 32bit, ~32 for 64bit, etc). This performs nomixing on the high (shift amount of) bits and mixes the high order bits with the loworder.
This means the last right xorshift of the mixing function only effects the low order bits of the result (the default of S2 is 33) so this can be dropped without effecting statistical quality if accounted for in using functions. Example the standard equispaced conversion to singles only uses the top 24 bits. Likewise integer results can be formed using only the high bits of the result.
The mixing functions constants are defined as: PRNS_MIX_S0
, PRNS_MIX_S1
, PRNS_MIX_S2
, PRNS_MIX_M0
, PRNS_MIX_M1
.
The final rightxorshift is disabled by defining: PRNS_NO_FINAL_XORSHIFT
The mixing function of the main sequence is additionally never called directly it is instead a macro expansion which defaults to: #define PRNS_MIX(X) prns_mix(X)
The down sequence likewise uses: #define PRNS_MIX_D(X) prns_mix(X)
The default provided bitmixing function is intended to be a middleground solution for people that do not want to bother with providing one. This default mixing functions performs respectably in the generalized test battery Crush^{3} (report) which is massive overkill.
If you usage requires less than a couple of hundred thousand random numbers per simulation step per entity then a lighter weight mixing function targeting Smallcrush^{4} is more appropriate: #define PRNS_SMALLCRUSH
sets the mixing function to a single right xorshiftmultiply (report).
These tests were run using the provided driver (SEE: Resources)
An alternate solution for (M)LCGs and PCGs
If we reduce the problem statement to focus on the ability to query the current member and either increment or decrement the position, then we could instead use a poweroftwo LCG. Then the forward state update looks like:
\[\begin{equation} \label{eq:lcgf} u_{i+1} = \left( a\thinspace u_{i} + b \right) \bmod 2^{w} \end{equation}\]reversing the recurrence gives us (dropping the explict mod):
\[\begin{equation} \label{eq:lcgr} u_{i1} = \left( a^{1}\left( u_{i}  b \right) \right) = \left( a^{1} u_{i}  a^{1}b \right) = \left( a_r\thinspace u_{i} + b_r \right) \end{equation}\]where $ a^{1} $ is the modulo inverse of $a$. So we have two standard poweroftwo LCG state updates where one simply walks the sequence in the opposite direction.
For the mixing function we could turn to one of those provide by O’Neill in her PCG paper^{5}. In this case would would have identical performance and statistical quality as the equivalent PCG with the added ability to walk both directions at the same cost.
The recurrence relation $ \eqref{eq:lcgf} $ can be expressed in closed form (without the modulo) for ‘k’ steps forward as:
\[\begin{equation} \label{eq:lcgcf} u_{i+k} = a^k\thinspace u_i + \left( \frac{a^k1}{a1} \right) b \end{equation}\]Equation $ \eqref{eq:lcgcf} $ can be evaluated using divideandconquer. If ‘k’ is negative we could still evaulate directly with the forward form but that would insure that we would be performing the maximum number of iterations. Effectively we would be walking the long way around the circle between the two points. Instead it is possible to use the reversed sequence constants. Obviously we can also compute the member at a given position using $ \eqref{eq:lcgcf} $ as well.
Note that the C++ version of PCG include skipping features, although it does not appear that the code takes advantage of being able to walk the minimal path.
Equation $ \eqref{eq:lcgcf} $ can be rewritten as:
\[u_{i+k} = M \thinspace u_i + C \thinspace b\]Where $M = a^k$ which has various well know divideandconquer forms and:
\[C = \left( \frac{a^k1}{a1} \right) = 1 + a + a^2 + ... + a^{k1}\]The basics are pretty straight forward and are implemented in: (lcgs.h)
This is currently defaulting to the dumbest mixing on the planet. I need to sit down an think about it. The source includes PCG mixing functions which you have to set macros to get included…not public domain.
Fabian Giesen has written a better description of sequence skipping for LCGs than this HERE with various choices of improved exp performance.
Generalization to any bijection state update
I would argue that all interesting state updates are bijections which can be expressed as the matrix equation over some field:
\[u_{i+1} = M\ u_{i}\]so $ M $ is invertible and there exists some scheme to produce $ M^{k} $, say by precomputing code for various $ k $ values. So there are O(1) formulations for reversed sequence walking and O(ln n) for skipping. However there are not any obvious choices such that forward/reverse walking are roughly the same complexity. As examples consider binarytogray code conversion vs. graytobinary or this paper^{6}.
Resources
References and Footnotes

Weyl sequence overview (local post) ↩

Better Bit Mixing  Improving on MurmurHash3’s 64bit Finalizer, David Stafford (blog post) ↩

Crush is a TestU01^{7} battery, roughly speaking passing this means you’re good for ~$2 \times 10^6$ samples per problem. ↩

Smallcrush is a TestU01^{7} battery, roughly speaking passing this means you’re good for ~$2 \times 10^5$ samples per problem. ↩

PCG: A Family of Simple Fast SpaceEfficient Statistically Good Algorithms for Random Number Generation, Melissa E. O’Neill, (paper) ↩

Efficient Jump Ahead for $\mathbb{F_2}$Linear Random Number Generators, Haramoto, Matsumoto, Nishimura, Panneton, L’Ecuyer, 2008. (PDF) ↩

TestU01: A C library for empirical testing of random number generators, Pierre L’Ecuyer and Richard Simard, 2007. (original source/paper) (hacked source) ↩ ↩^{2}