- Generated by
1.8.13
|
LatMRG Guide
1.0
A software package to test and search for new linear congruential random number generators
|
This page covers the theory specific to LatMRG in two main sections and a third secondary one.
First, we present all the generators engines covered by the library. Along with the generators definitions, we present the equivalences that make it possible to study all those generatorse in LatMRG. In the second section, we present the construction of lattices for such random generators, as well as the computations to study such lattices. In the last section, we present some issues arrising when choosing parameters for a generator. This is an optionnal section that can be helpfull when starting searches with LatMRG.
We assume some familiarity with the concept of random number generator from from the reader. If this is not the case, we recommend reading the beginning of [rLEC17h]. Although this is not the shortest introduction to random number generators, it is very well suited for beginners and it covers both historical and modern issues of random number generation.
We first present the oldest and simplest generator engine that LatMRG covers. Given \(a,m,c \in \mathbb{Z}\), usually also with the constraints \(0 < a < m\) and \(0 \leq c < m,\ 1 \equiv c \pmod 2\), we can get pseudo-random number in \((0,1)\) by generating the sequence \(\{x_n\}\)
\[ x_n = a x_{n-1} + c \pmod m \]
and by taking the output of \(u_n = x_n / m\). Note that if \(c = 0\), the components also need to be such that \(\gcd(a,m) = 1\).
This kind of generator is called a linear congruential generator LCG. It is easy to implement and runs fast. But as computing power grew, It became clear that it was not practical to have a LCG as a random number generator. These generator have a high correlation between their points. It is also impractical to make on with a long enough period as the maximal possible period for a LCG is \(m-1\).
The natural extension of this engine that arose was the augmentation of the order of the recurrence as follows:
\[ x_n = a_1 x_{n-1} + \cdots + a_k x_{n-k} \pmod m \]
where \(a_k \neq 0\). This engine can output random numbers in the same way by taking \(u_n = x_n / m\). Using this kind of recurrence mostly solves the problems of LCGs as it increases the states space considerably. This kind of engine is called a multiple recursive generator (MRG). Note that LCG generators simply are a specific case of MRG generators. It turns out that most of the generators presented after this point are equivalent to a MRG in some way and can be considered simply as MRGs when developping the theory. This is for this reason that this software is called LatMRG. Lattices and Merit presents the main results to study a MRG, but they also have been broadly studied in the litterature in, for example, [rLEC93a] [rLEC88a], and [rLEC97b].
Using a MRG is much slower that using a LCG so many avenues have been studied as a way to generate random numbers on a period large enough to be useable, but that is faster than an MRG. This is what [rMAR91a] introduced. These engines are fast and also have long periods and use the following recurrences with \(r > s\) and \(b \in \mathbb{Z}_{>0}\):
\[ \begin{array}{llll} \text{(AWC)} & x_n & = & x_{n-s} + x_{n-r} + c_n \pmod b, \\ & c_{n+1} & = & \mathbb{1}(x_{n-s} + x_{i-r} + c_n \geq b). \\ \text{(AWC-c)} & x_n & = & -x_{n-s} - x_{n-r} - c_n - 1 \pmod b,\\ & c_{n+1} & = & \mathbb{1}(x_{n-s} + x_{i-r} + c_n \geq b).\\ \text{(SWB-I)} & x_n & = & x_{n-s} - x_{n-r} - c_n \pmod b,\\ & c_{n+1} & = & \mathbb{1}(x_{n-s} - x_{n-r} - c_n < 0).\\ \text{(SWB-II)} & x_n & = & - x_{n-s} + x_{n-r} - c_n \pmod b,\\ & c_{n+1} & = & \mathbb{1}(x_{n-s} - x_{n-r} - c_n < 0).\\ \end{array} \]
These generators are called Add-with-Carry (AWC) and Subtract-with-Borrow (SWB) (we follow the naming scheme of [rTEZ93a] for the variants). The fact that these generators do not need multiplication operations mean that they will be very fast with b as a power of 2. With a good choice of parameters, the period of the recurrence can be up to around \(b^r\), which means that this generator solves both the period length problem and the speed one.
The sequence \(\{x_n\}\) can be used to generate random \(\mathcal{U}(0,1)\) numbers with an arbitrary number of digits in base \(b\) with
\[ u_n = \sum_{j=1}^L x_{Ln + j + 1} b^{-j}. \]
By taking \(M\) as follows,
| Generator | M |
|---|---|
| AWC | \(b^r+b^s-1\) |
| AWC-c | \(b^r+b^s+1\) |
| SWB-I | \(b^r-b^s+1\) |
| SWB-II | \(b^r-b^s-1\) |
We can build an equivalence between AWC/SWB generators and the following LCG ([rTEZ93a])
\begin{align} X_n & = AX_{n-1} \pmod M \\ v_n & = X_n / M \\ w_n & = X_{Ln}/M. \end{align}
When we say equivalence, this is in the sense that, with the right initial state, the output \(w_i\) and \(u_i\) have the same decimal digits in base \(b\) up to the precision of \(u_i\). That is: \(u_i = b^{-L} \lfloor b^L w_i \rfloor\).
It turns out that AWC/SWB generators have a few faults in theoretical tests. [rMAR94a] suggested a similar kind of recurrence to solve this issue called multiply-with-carry (MWC)
\[ \begin{array}{lll} x_n & = & (a_1 x_{n-1} + \cdots + a_k x_{n-k} + c_{n-1})d \pmod b,\\ c_n & = & \lfloor (a_0 x_n + a_1 x_{n-1} + \cdots + a_k x_{n-k})/b\rfloor, \end{array} \]
where \(\gcd(a_0, b) = 1\) and \(1 \equiv -a_0d \pmod b\). This engine uses the same output function as AWC/SWB generators \( u_n = \sum_{j=1}^L x_{Ln + j + 1} b^{-j} \) and the results of [rTEZ93a] can be tweaked as in [rCOU97a]. Take \(M = \sum_{l = 0}^k a_k b^l\) and \(0 < A < m\) with \(1 \equiv Ab \pmod m\), then the same recurrence
\begin{align} X_n & = AX_{n-1} \pmod M \\ v_n & = X_n / M \\ w_n & = X_{Ln}/M. \end{align}
is equivalent to the MWC generator with the same precision: \(u_i = b^{-L} \lfloor b^L w_i \rfloor\).
Generators using a carry can easily be constructed to have a large period and to operate very fast, but the fact that they are equivalent to LCGs means that their structure is most of the time flawed, especially in 2 dimensions. To solve this problem, it was suggested to use multiple MRGs and to combine their output. This reduces the computing cost of big MRGs, while obtaining much better distribution than carry generators.
Let sequences \(\{\{x_n^{(1)}\},\{x_n^{(2)}\},\ldots,\{x_n^{(\ell)}\}\}\) and \(\{\{u_n^{(1)}\},\{u_n^{(2)}\},\ldots,\{u_n^{(\ell)}\}\}\) be produced by MRGs with \(a_{i,j}\) the j-th coefficient of the i-th recurrence, \(k_i\) the order of the i-th recurrence and \(m_i\) the modulus of the i-th recurrence. It is possible to get pseudo-random numbers as:
\[ \begin{array}{lll} \tilde{u}_n & = & \left(\sum_{i = 1}^\ell \delta_{i} x_n^{(i)} \pmod {m_1} \right) / m_1 \\ w_n & = & \sum_{i = 1}^\ell \frac{\delta_{i} x^{(i)}_n}{m_i} \pmod 1. \end{array} \]
[rLEC96b] presents a way to build a MRG that has an output \(\{u_n\}\) identical to \(\{w_n\}\) if their first \(k = \max(k_1, \ldots, k_\ell)\) output values are the same. Define
\begin{align} m & = \prod_{i = 1}^\ell m_i \\ b & = \sum_{i=1}^\ell \frac{\delta_i b_i m}{m_i} \pmod m \\ n_j & \text{with } 1 \equiv n_j (m /m_j) \pmod {m_j},\ 1 \leq j \leq \ell \\ a_j & = \sum_{i=1}^\ell \frac{a_{i,j} n_i m}{m_i} \pmod m,\ 1 \leq j \leq k. \end{align}
Then if \((u_0, \ldots, u_{k-1}) = (w_0, \ldots, w_{k-1})\), then \(w_n = u_n\) for all \(n\in\mathbb{N}\) given
\begin{align} x_n & = a_1 x_{n-1} + \cdots + a_k x_{n-k} + b \pmod m \\ u_n & = x_n/m. \end{align}
The final type of generator that we present is slightly different. This is because it cannot be transformed to an equivalent MRG. Instead, this is a generalization of an MRG.
Take \(\mathbf{A}\in \mathbb{Z}^{k\times k}\) and \(s_n \in \mathbb{Z}^k\). We can generate a random number vector with
\[ s_n = \mathbf{A} s_{n-1} \pmod m \]
and \(\mathbf{u}_n = (1/m) \cdot s_n\). This is what we call a matrix multiple recursive generator (MMRG). This kind of generator can also be used to generate any size of vector, including single random numbers. Let \(\mathbf{u}_n(i)\) be the i-th component for vector \(\mathbf{u}_n\) and \(p = qk + r,\ 0 \leq r < k\) be the size of vectors wanted. We can get the random vector sequence \(\{v_n = (\mathbf{u}_{n(q+1)}(1), \ldots, \mathbf{u}_{n(q+1)}(k), \ldots, \mathbf{u}_{n(q+1)+q}(1), \ldots, \mathbf{u}_{n(q+1)+q}(r)) \in \mathbb{Z}^p\}\).
We will see below that a MRG of order \(k\) with multipliers \((a_1, \ldots, a_k)\) and modulo \(m\) has the same output structure as a MMRG with
\[ \mathbf{A} = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \ddots & 0 \\ \vdots & & \ddots & \ddots & \vdots \\ 0 & \cdots & \cdots & 0 & 1 \\ a_1 & a_2 & \cdots & \cdots & a_k \end{bmatrix}^k \]
As explained in Theory on Lattices, lattices of dimension \(t\) are discrete subgroups of \(\mathbb{R}^t\) generated by the integer combinations of a set of linearly independant (over the integers) vectors. We now explain how random number generators of the MRG familly can be studied from the lattice of their points, and how to obtain the lattice of these generators.
Take
\[\Psi_t = \{ \mathbf{u}_0 =(x_0 /m, \ldots, x_{t-1}/m): (x_0, \dots, x_{k-1}) \in \mathbb{Z}_m^k\}\]
for a MRG of order \(k\). This is the set of all first \(t\) components vectors the generator will output. The periodic contination of this set is a lattice. We will name this lattice \(L_t\), thus
\[ \Psi_t + \mathbb{Z}^t = L_t. \]
This lattice has basis
\[ \mathbf{V} = \begin{pmatrix} 1/m & 0 & \dots & 0 & x_{1,k+1} & \dots & x_{1,t} \\ 0 & 1/m & \dots & 0 & x_{2,k+1} & \dots & x_{2,t} \\ \vdots & \vdots & \ddots & \vdots & \vdots & & \vdots \\ 0 & 0 & \dots & 1/m & x_{k,k+1} & \dots & x_{k,t} \\ 0 & 0 & \dots & 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 & 0 & \dots & 1 \end{pmatrix} \]
and its dual has basis
\[ \mathbf{W} = \begin{pmatrix} m & 0 & \dots & 0 & 0 & \dots & 0 \\ 0 & m & \dots & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & & \vdots \\ 0 & 0 & \dots & m & 0 & \dots & 0 \\ -x_{1,k+1} & -x_{2,k+1} & \dots & -x_{k,k+1} & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ -x_{1,t} & -x_{2,t} & \dots & -x_{k,t} & 0 & \dots & 1 \end{pmatrix}. \]
Where \(x_{i,j}\) is the \(j\)-th x value of the generator with initial state \( \mathbf{e}_i = (0, \dots, 1, \dots, 0) \in \mathbb{Z}^k\), the \(i\)-th canonical vector in dimension \(k\). LatMRG can build this lattice for any of the generators mentioned.
Instead of forming vectors with successive values like in the above definition of \(\Psi_t\), one can form vectors with values that are some distance apart in the sequence (so-called "leapfrog" values). a set of fixed integers. Define
\[ \Psi_t(I) = \{(u_{i_1},…,u_{i_t}) \mid(x_0,…,x_{k-1}) \in\mathbb{Z}_m^k\} \]
and let \(L_t(I) = \Psi_t(I) + \mathbb{Z}^t\). If we assume that \(0\le i_1 < i_2 < \cdots< i_t\), this \(L_t(I)\) is the projection of the lattice \(L_{i_t+1}\) over the \(t\)-dimensional subspace determined by the coordinates that belong to \(I\). For \((i_1,…,i_t) = (0,…,t-1)\), one has \(L_t(I) = L_t\). Using what is available in LatticeTester one can build a basis for \(L_t\) and its dual in this more general case. This means LatMRG can perform the usual lattice analysis on these kind of lattice also. Further details and examples are given in [rLEC97c].
The lattice structure also means that all points of \(L_t\) lie in a family of equidistant parallel hyperplanes. Examining geometrical properties of these hyperplanes provides information on the uniformity of of the random numbers generated. Another important property that can be used to study of MRGs is their period length. Although a long period does not garantee good uniformity, it is ideal for better distribution properties.
Among all such families hyperplanes that cover all the points, choose the one for which the successive hyperplanes are farthest apart. The distance between these successive hyperplanes is equal to \(1/\ell_t\) where \(\ell_t\) is the Euclidean length of the shortest nonzero vector in the dual lattice \(L_t^*\). It is possible to normalize this value, such as to obtain the spectral test for the lattice.
To build figures of merit from the spectral test, we take the worst-case value of \(\ell_t/\ell_t^*(m^k)\) over certain values of \(t\) and for selected projections on lower-dimensional subspaces. More specifically, let \(\ell_I\) denote the length of the shortest nonzero vector \(\mathbf{v}\) in \(L_t^*(I)\), and \(\ell_t = \ell_{\{1,…,t\}}\) as before. For arbitrary positive integers \(t_1\ge\cdots\ge t_d \ge d\), consider the worst-case figure of merit
\[ M_{t_1,…,t_d} = \min\left[ \min_{k+1\le t\le t_1} \ell_t/\ell_t^*(m^k),\min_{2\le s\le k} \min_{I\in S(s,t_s)} \ell_I/m, \min_{k+1\le s\le d} \min_{I\in S(s,t_s)} \ell_I/\ell_s^*(m^k) \right], \]
where \(S(s,t_s) = \{I=\{i_1,…,i_s\} \mid1 = i_1 < \cdots< i_s\le t_s\}\). This figure of merit makes sure that the lattice is good in projections over \(t\) successive dimensions for all \(t\le t_1\), and over non-successive dimensions that are not too far apart. Note that when \(s\le k\), the smallest distance between hyperplanes that can be achieved in \(s\) dimensions for the MRG is \(1/m\) because it is possible for the generator to span all points with coordinates \(z/m\) with \(z\in \mathbb{Z}\).
The figure of merit \(M_{t_1} = \min_{2\le s\le t_1} \ell_s/\ell_s^*(n)\) (with \(d=1\)) has been widely used for ranking and selecting LCGs and MRGs [sFIS96a], [rLEC99b], [rLEC99c]. The quantity \(M_{t_1,…,t_d}\) is a worst case over \((t_1-d) + \sum_{s=2}^d \binom{t_s-1}{s-1}\) projections, and this number increases quickly with \(d\) unless the \(t_s\) are very small.
When too many projections are considered, there are inevitably some that are bad. The worst-case figure of merit is (practically) always small, and cannot distinguish between good and mediocre behavior in the most important projections. Moreover, the time to compute \(M_{t_1,…,t_d}\) increases with the number of projections. We therefore suggest considering only the projections deemed important when first comparing generators. It is then possible to compare good generators for those projections more extensively. We suggest using the criterion with \(d\) equal to 4 or 5, and \(t_s\) decreasing with \(s\).
Instead of considering the shortest nonzero vector in the dual lattice, one can consider the shortest nonzero vector in the primal lattice \(L_t\). Its length represents the distance to the nearest other lattice point from any point of the lattice. A small value means that many points are placed on the same line, at some fixed distance apart.
Another way of measuring the quality of a lattice is in terms of the relative lengths of the smallest and largest vectors in a reduced basis. It has been historicaly important to use this ratio for a Minkowski-reduced basis (MRLB) (see [rAFF85a], [rAFF88a], [rGRO88a] for more details). Roughly, a MRLB is a basis for which the vectors are in some sense the most orthogonal.
The ratio of the sizes of the shortest and longest vectors of a MRLB is called its Beyer-quotient. In general, a given lattice may have several MRLBs, all with the same length of the shortest vector, but perhaps with different lengths of the longest vector, and thus different Beyer quotients. We define \(q_t(I)\) as the maximum of the Beyer quotients of all MRLBs of \(L_t(I)\), and denote \(q_t(\{1,…,t\})\) by \(q_t\). We prefer values of \(q_t(I)\) close to 1. Similar to \(M_{t_1,…,t_d}\), we define
\[ Q_{t_1,…,t_d} = \min\left[ \min_{k+1\le t\le t_1} q_t,\min_{2\le s\le d} \min_{I\in S(s,t_s)} q_t(I) \right]. \tag{Q} \]
Since computing \(q_t\) is much more time consuming than computing the spectral test and this mesure is not unique, this test is seldom used nowadays.
The choice of modulo is a recurrent source of discussion when building such generators. This is because a modulo operation takes a variable time for variable modulus. Most of the time, there is a compromise to do between simple modulus like \(2^q\) (which can be computed with a simple xor operation) and \(p\) big prime number. The former is much faster to compute, but the latter usually performs much better in theoretical and statistical tests.
LatMRG can deal with very large moduli and multipliers. There is no limit on size other than the size of the computer memory (and the CPU time). For example, a generator with a modulus of a few hundred bits can be analyzed easily. Operations on large integers are performed using the GNU multi-precision package GMP [iGMP06a]. GMP is a portable library written in C for arbitrary precision arithmetic on integers, rational numbers, and floating-point numbers. For vectors, matrices of large numbers and polynomials, we use NTL [iSHO05a]. NTL is a high-performance, portable C++ library providing data structures and algorithms for manipulating arbitrary length integers, and for vectors, matrices, and polynomials over the integers and over finite fields. NTL uses GMP as an underlying package for dealing with large numbers.
Of course, arithmetic operations with these structures are performed in software and are significantly slower than the standard operations supported by hardware. For this reason, most of the basic (low-level) operations required by our higher-level classes have been implemented in two or three versions. When building a basis or checking maximal period conditions, the modulus and multipliers can be represented either as int64_t’s (64-bit integers) or ZZ’s (arbitrary large integers). After a lattice basis and its dual have been constructed, when working on the basis (finding a shortest vector, Minkowski reduction, etc.), the vector elements can be represented either as double’s (64-bit floating-point numbers) or RR’s (arbitrary large floating-point numbers).
When performing a search for good generators, for instance, one can first perform all the "screening" computations (involving many generators) using standard type int64_t, and then recompute (verify) with the large floating-point numbers RR only for the retained generators.
Beware! A lot of the manipulations needed for 64 bits generators need more that 64 bits. Using int64_t to search and test 64 bits generators might cause overflows.
1.8.13