Subsections

• 2.4.1 A General Framework
• 2.4.2 Measures of Uniformity
• 2.4.3 Lattice Structure in Spaces of Polynomials and Formal Series
• 2.4.4 The LFSR Generator
• 2.4.5 The GFSR and Twisted GFSR
• 2.4.6 Combined Linear Generators Over $\mathbb{F}_2$

2.4 Generators Based on Recurrences Modulo 2

2.4.1 A General Framework

It seems natural to exploit the fact that computers work in binary arithmetic and to design RNGs defined directly in terms of bit strings and sequences. We do this under the following framework, taken from [54]. Let ${\mathbb{F}}_2$ denote the finite field with two elements, 0 and $1$, in which the operations are equivalent to addition and multiplication modulo $2$. Consider the RNG defined by a matrix linear recurrence over  ${\mathbb{F}}_2$, as follows:

$$\boldsymbol{x}_i = \boldsymbol{A} \boldsymbol{x}_{i-1}\;,$$ (2.12)

$$\boldsymbol{y}_i = \boldsymbol{B} \boldsymbol{x}_{i}\;,$$ (2.13)

$$u_i = \sum_{\ell=1}^w y_{i,\ell-1} 2^{-\ell} = y_{i,0}\; y_{i,1}\; y_{i,2}\; \cdots\;,$$ (2.14)

where $\boldsymbol{x}_i = (x_{i,0},\ldots,x_{i,k-1})^{{{\sf T}}} \in {\mathbb{F}}_2^k$ is the $k$-bit state vector at step $i$, $\boldsymbol{y}_i = (y_{i,0},\ldots,y_{i,w-1})^{{{\sf T}}} \in {\mathbb{F}}_2^w$ is the $w$-bit output vector at step $i$, $k$ and $w$ are positive integers, $\boldsymbol{A}$ is a $k\times k$ transition matrix with elements in ${\mathbb{F}}_2$, $\boldsymbol{B}$ is a  $w \times k$ output transformation matrix with elements in ${\mathbb{F}}_2$, and $u_i\in [0,1)$ is the output at step $i$. All operations in (2.12) and (2.13) are performed in  ${\mathbb{F}}_2$.

It is well-known ([80,36]) that when the $\boldsymbol{x}_i$'s obey (2.12), for each $j$, the sequence $\{x_{i,j},\, i\ge 0\}$ follows the linear recurrence

$$x_{i,j} = \left(\alpha_1 x_{i-1,j} + \cdots + \alpha_k x_{i-k,j}\right) \mod 2\;,$$ (2.15)

whose characteristic polynomial $P(z)$ is the characteristic polynomial of  $\boldsymbol{A}$, i.e.,

$$P(z) = (\boldsymbol{A}-z\boldsymbol{I}) = z^k - \alpha_1 z^{k-1} - \cdots - \alpha_{k-1} z - \alpha_k\;,$$

where $\boldsymbol{I}$ is the identity matrix and each $\alpha_j$ is in ${\mathbb{F}}_2$. The sequences $\{y_{i,j},\, i\ge 0\}$, for $0\le j < w$, also obey the same recurrence (although some of them may follow recurrences of shorter order as well in certain situations, depending on $\boldsymbol{B}$). We assume that $\alpha_k=1$, so that the recurrence (2.15) has order $k$ and is purely periodic. Its period length is $2^k-1$ (i.e., maximal) if and only if $P(z)$ is a primitive polynomial over ${\mathbb{F}}_2$ ([80,31]).

To jump ahead directly from $\boldsymbol{x}_i$ to $\boldsymbol{x}_{i+\nu}$ with this type of generator, it suffices to precompute the matrix $\boldsymbol{A}^\nu$ (in ${\mathbb{F}}_2$) and then multiply $\boldsymbol{x}_i$ by this matrix.

Several popular classes of RNGs fit this framework as special cases, by appropriate choices of the matrices $\boldsymbol{A}$ and  $\boldsymbol{B}$. This includes the Tausworthe or LFSR, polynomial LCG, GFSR, twisted GFSR, Mersenne twister, multiple recursive matrix generators, and combinations of these ([54,78,81,88]). We detail some of them after discussing measures of their uniformity.

2.4.2 Measures of Uniformity

The uniformity of point sets $\Psi_I$ produced by RNGs based on linear recurrences over ${\mathbb{F}}_2$ is usually assessed by measures of equidistribution defined as follows ([38,54,45,88]). For an arbitrary vector $\boldsymbol{q} = (q_1,\ldots,q_t)$ of non-negative integers, partition the unit hypercube $[0,1)^t$ into $2^{q_j}$ intervals of the same length along axis $j$, for each $j$. This determines a partition of $[0,1)^t$ into $2^{q_1 + \cdots + q_t}$ rectangular boxes of the same size and shape. We call this partition the $\boldsymbol{q}$-equidissection of the unit hypercube.

For some index set $I = \{i_1,\ldots,i_t\}$, if $\Psi_I$ has $2^k$ points, we say that $\Psi_I$ is $\boldsymbol{q}$-equidistributed in base $2$ if there are exactly $2^{q}$ points in each box of the $\boldsymbol{q}$-equidissection, where $k-q = q_1 + \cdots + q_t$. This means that among the $2^k$ points $(x_{j_1},\ldots,x_{j_t})$ of $\Psi_I$, if we consider the first $q_1$ bits of $x_{j_1}$, the first $q_2$ bits of $x_{j_2}, \ldots$, and the first $q_t$ bits of $x_{j_t}$, each of the $2^{k-q}$ possibilities occurs exactly the same number of times. This is possible only if $q\le k$.

The $\boldsymbol{q}$-equidistribution of $\Psi_I$ depends only on the first $q_j$ bits of $x_{i_j}$ for $1\le j\le t$, for the points $(x_{i_1},\ldots,x_{i_t})$ that belong to $\Psi_I$. The vector of these $q_1 + \cdots + q_t = k-q$ bits can always be expressed as a linear function of the $k$ bits of the initial state $\boldsymbol{x}_0$, i.e., as $\boldsymbol{M}_{\boldsymbol{q}} \boldsymbol{x}_0$ for some $(k-q)\times k$ binary matrix $\boldsymbol{M}_{\boldsymbol{q}}$, and it is easily seen that $\Psi_I$ is $\boldsymbol{q}$-equidistributed if and only if $\boldsymbol{M}_{\boldsymbol{q}}$ has full rank $k-q$. This provides an easy way of checking equidistribution ([22,38,88]).

If $\Psi_I$ is $(\ell,\ldots,\ell)$-equidistributed for some $\ell\ge 1$, it is called $t$-distributed with $\ell$ bits of accuracy, or $(t,\ell)$-equidistributed ([38]). The largest value of $\ell$ for which this holds is called the resolution of the set $\Psi_I$ and is denoted by $\ell_I$. This value has the upper bound $\ell_t^{\ast} = \min(\lfloor k/t\rfloor,\, w)$. The resolution gap of $\Psi_I$ is defined as $\delta_I = \ell_t^{\ast} - \ell_I$. In the same vein as for MRGs, a worst-case figure of merit can be defined here by

$$\Delta_{\mathcal{J}} = \max_{I\in{\mathcal{J}}} \delta_I\;,$$

where ${\mathcal{J}}$ is a preselected class of index sets $I$.

The point set $\Psi_I$ is a $(q,k,t)$-net in base $2$ (often called a $(t,m,s)$-net in the context of quasi-Monte Carlo methods, where a different notation is used ([80])), if it is $(q_1,\ldots,q_t)$-equidistributed in base $2$ for all non-negative integers $q_1,\ldots,q_t$ summing to $k-q$. We call the smallest such $q$ the $q$-value of $\Psi_I$. The smaller it is, the better. One candidate for a figure of merit could be the $q$-value of $\Psi_t$ for some large $t$. A major drawback of this measure is that it is extremely difficult to compute for good long-period generators (for which $k-q$ is large), because there are too many vectors $\boldsymbol{q}$ for which equidistribution needs to be checked. In practice, one must settle for figures of merit that involve a smaller number of equidissections.

When $\delta_I = 0$ for all sets $I$ of the form $I = \{0,\ldots,t-1\}$, for $1\le t\le k$, the RNG is said to be maximally equidistributed or asymptotically random for the word size $w$ ([38,88,91]). This property ensures perfect equidistribution of all sets $\Psi_t$, for any partition of the unit hypercube into subcubes of equal sizes, as long as $\ell\le w$ and the number of subcubes does not exceed the number of points in $\Psi_t$. Large-period maximally equidistributed generators, together with their implementations, can be found in [43,54], and [84], for example.

2.4.3 Lattice Structure in Spaces of Polynomials and Formal Series

The RNGs defined via (2.12)-(2.14) do not have a lattice structure in the real space like MRGs, but they do have a lattice structure in a space of formal series, as explained in [12,45,65], and [88]. The real space ${\mathbb{R}}$ is replaced by the space ${\mathbb{L}}_2$ of formal power series with coefficients in ${\mathbb{F}}_2$, of the form $\sum_{\ell=\omega}^\infty x_\ell z^{-\ell}$ for some integer $\omega$. In that setting, the lattices have the form

$${\mathcal{L}}_t = \left\{\boldsymbol{v}(z) = \sum_{j=1}^t h_j(z) ... ...v}_j(z)\quad \text{such that each}\quad h_j(z) \in {\mathbb{F}}_2[z]\right\}\;,$$

where ${\mathbb{F}}_2[z]$ is the ring of polynomials with coefficients in ${\mathbb{F}}_2$, and the basis vectors  $\boldsymbol{v}_j(z)$ are in ${\mathbb{L}}_2^t$. The elements of the dual lattice ${\mathcal{L}}_t^{\ast}$ are the vectors $\boldsymbol{h}(z)$ in ${\mathbb{L}}_2^t$ whose scalar product with any vector of ${\mathcal{L}}_t$ is a polynomial (in ${\mathbb{F}}_2[z]$). We define the mapping $\varphi : {\mathbb{L}}_2\to{\mathbb{R}}$ by

$$\varphi\left(\sum_{\ell=\omega}^\infty x_\ell z^{-\ell}\right) = \sum_{\ell=\omega}^\infty x_\ell 2^{-\ell}\;.$$

Then, it turns out that the point set $\Psi_t$ produced by the generator is equal to $\varphi({\mathcal{L}}_t)\cap [0,1)^t$. Moreover, the equidistribution properties examined in Sect. 2.4.2 can be expressed in terms of lengths of shortest vectors in the dual lattice, with appropriate definitions of the length (or norm). Much of the theory and algorithms developed for lattices in the real space can be adapted to these new types of lattices ([12,45,65,88]).

2.4.4 The LFSR Generator

The Tausworthe or linear feedback shift register (LFSR) generator ([87,38,88]) is a special case of (2.12)-(2.14) with $\boldsymbol{A} = \boldsymbol{A}_0^s$ (in ${\mathbb{F}}_2$) for some positive integer $s$, where

$$\boldsymbol{A}_0 = \left( \begin{tabular}{cccc} & $1$\ & & \\ & &... ...& & & $1$\ \\ $a_k$\; &$a_{k-1}$\ & $\ldots$\ & $a_1$\ \end{tabular} \right)\;,$$ (2.16)

$a_1,\ldots,a_k$ are in ${\mathbb{F}}_2$, $a_k=1$, and all blank entries in the matrix are zeros. We take $w \le k$ and the matrix $\boldsymbol{B}$ contains the first $w$ lines of the $k\times k$ identity matrix. The RNG thus obtained can be defined equivalently by

$$x_i = a_1 x_{i-1} + \cdots + a_k x_{i-k} \mod 2\;,$$ (2.17)

$$u_i = \sum_{\ell=1}^w x_{is+\ell-1} 2^{-\ell}\;.$$ (2.18)

Here, $P(z)$ is the characteristic polynomial of the matrix $\boldsymbol{A}_0^s$, not the characteristic polynomial of the recurrence (2.17), and the choice of $s$ is important for determining the quality of the generator. A frequently encountered case is when a single $a_j$ is nonzero in addition to $a_k$; then, $P(z)$ is a trinomial and we say that we have a trinomial-based LFSR generator. These generators are known to have important statistical deficiencies ([77,88]) but they can be used as components of combined RNGs.

LFSR generators can be expressed as LCGs in a space of polynomials ([89,88,36]). With this representation, their lattice structure as discussed in Sect. 2.4.3 follows immediately.

2.4.5 The GFSR and Twisted GFSR

Here we take $\boldsymbol{A}$ as the $pq \times pq$ matrix

$$\boldsymbol{A} = \left( \begin{tabular}{ccccc} & & $ \boldsymbol{... ... & \\ & & $\ddots$\ & & \\ & & & $\boldsymbol{I}_p$\ & \\ \end{tabular} \right)$$

for some positive integers $p$ and $q$, where $\boldsymbol{I}_p$ is the $p \times p$ identity matrix, $\boldsymbol{S}$ is a $p \times p$ matrix, and the matrix $\boldsymbol{I}_p$ on the first line is in columns $(r-1)p+1$ to $rp$ for some positive integer $r$. Often, $w=p$ and $\boldsymbol{B}$ contains the first $w$ lines of the $pq \times pq$ identity matrix. If $\boldsymbol{S}$ is also the identity matrix, the generator thus obtained is the trinomial-based generalized feedback shift register (GFSR), for which $\boldsymbol{x}_i$ is obtained by a bitwise exclusive-or of $\boldsymbol{x}_{i-r}$ and $\boldsymbol{x}_{i-q}$. This gives a very fast RNG, but its period length cannot exceed $2^q-1$, because each bit of $\boldsymbol{x}_i$ follows the same binary recurrence of order $k=q$, with characteristic polynomial $P(z) = z^q - z^{q-r} - 1$.

More generally, we can define $\boldsymbol{x}_i$ as the bitwise exclusive-or of $\boldsymbol{x}_{i-r_1}, \boldsymbol{x}_{i-r_2}, \ldots, \boldsymbol{x}_{i-r_d}$ where $r_d=q$, so that each bit of $\boldsymbol{x}_i$ follows a recurrence in ${\mathbb{F}}_2$ whose characteristic polynomial $P(z)$ has $d+1$ nonzero terms. However, the period length is still bounded by $2^q-1$, whereas considering the $pq$-bit state, we should rather expect a period length close to $2^{pq}$. This was the main motivation for the twisted GFSR (TGFSR) generator. In the original version introduced by [75], $w=p$ and the matrix $\boldsymbol{S}$ is defined as the transpose of $\boldsymbol{A}_0$ in (2.16), with $k$ replaced by $p$. The characteristic polynomial of $\boldsymbol{A}$ is then $P(z) = P_S(z^q+z^m)$, where $P_S(z) = z^p - a_p z^{p-1} - \cdots - a_1$ is the characteristic polynomial of $S$, and its degree is $k = pq$. If the parameters are selected so that $P(z)$ is primitive over  ${\mathbb{F}}_2$, then the TGFSR has period length $2^k-1$. [76] pointed out important weaknesses of the original TGFSR and proposed an improved version that uses a well-chosen matrix  $\boldsymbol{B}$ whose lines differ from those of the identity. The operations implemented by this matrix are called tempering and their purpose is to improve the uniformity of the points produced by the RNG. The Mersenne twister ([78,83]) is a variant of the TGFSR where $k$ is slightly less than $pq$ and can be a prime number. A specific instance proposed by [78] is fast, robust, has the huge period length of $2^{19937}-1$, and has become quite popular.

In the multiple recursive matrix method of [81], the first row of $p \times p$ matrices in $\boldsymbol{A}$ contains arbitrary matrices. However, a fast implementation is possible only when these matrices are sparse and have a special structure.

2.4.6 Combined Linear Generators Over $\mathbb{F}_2$

Many of the best generators based on linear recurrences over ${\mathbb{F}}_2$ are constructed by combining the outputs of two or more RNGs having a simple structure. The idea is the same as for MRGs: select simple components that can run fast but such that their combination has a more complicated structure and highly-uniform sets $\Psi_I$ for the sets $I$ considered important.

Consider $J$ distinct recurrences of the form (2.12)-(2.13), where the $j$th recurrence has parameters $(k, w, \boldsymbol{A}, \boldsymbol{B}) = (k_j, w, \boldsymbol{A}_j, \boldsymbol{B}_j)$ and state $\boldsymbol{x}_{j,i}$ at step $i$, for $j=1,\ldots,J$. The output of the combined generator at step $i$ is defined by

$$\boldsymbol{y}_i = \boldsymbol{B}_1\boldsymbol{x}_{1,i} \oplus \cdots \oplus\boldsymbol{B}_J\boldsymbol{x}_{J,i}\;,$$

$$u_i = \sum_{\ell=1}^w y_{i,\ell-1} 2^{-\ell}\;,$$

where $\oplus$ denotes the bitwise exclusive-or operation. One can show ([88]) that the period length $\rho$ of this combined generator is the least common multiple of the period lengths $\rho_j$ of its components. Moreover, this combined generator is equivalent to the generator (2.12)-(2.14) with $k=k_1 + \cdots + k_J$, $\boldsymbol{A} =$ diag $(\boldsymbol{A}_1,\ldots,\boldsymbol{A}_J)$, and $\boldsymbol{B} = (\boldsymbol{B}_1,\ldots, \boldsymbol{B}_J)$.

With this method, by selecting the parameters carefully, the combination of LFSR generators with characteristic polynomials $P_1(z),\ldots,P_J(z)$ gives yet another LFSR with characteristic polynomial $P(z) = P_1(z)\cdots P_J(z)$ and period length equal to the product of the period lengths of the components ([89,97,38,88]). Tables and fast implementations of maximally equidistributed combined LFSR generators are given in [38].

The TGFSR and Mersenne twister generators proposed in [76,78] and [83] cannot be maximally equidistributed. However, concrete examples of maximally equidistributed combined TGFSR generators with period lengths near $2^{466}$ and $2^{{1250}}$ are given in [54]. These generators have the additional property that the resolution gaps $\delta_I$ are zero for a class of small sets $I$ with indices not too far apart.