Improved Long-Period Generators Based on
Linear Recurrences Modulo 2FRANÇOIS PANNETON and PIERRE L’ECUYER
Université de Montréaland
MAKOTO MATSUMOTO
Hiroshima University

Fast uniform random number generators with extremely long periods have been defined and implemented based on linear recurrences modulo 2. The twisted GFSR and the Mersenne twister are
famous recent examples. Besides the period length, the statistical quality of these generators is
usually assessed via their equidistribution properties. The huge-period generators proposed so far
are not quite optimal in this respect. In this article, we propose new generators of that form with better equidistribution and “bit-mixing” properties for equivalent period length and speed. The stateof our new generators evolves in a more chaotic way than for the Mersenne twister. We illustrate
how this can reduce the impact of persistent dependencies among successive output values, which
can be observed in certain parts of the period of gigantic generators such as the Mersenne twister.
Categories and Subject Descriptors: G.4 [Mathematics of Computing]: Mathematical Software—
Algorithm design and analysis; I.6 [Computing Methodologies]: Simulation and Modeling
General Terms: Algorithms, Performance
Additional Key Words and Phrases: Random number generation, linear feedback shift register,
GFSR linear recurrence modulo 2, Mersenne twister

1. INTRODUCTION
Most uniform random number generators (RNGs) used in computational statistics and simulation are based on linear recurrences modulo 2, or modulo a large
This work has been supported by the Natural Sciences and Engineering Research Council of Canada(NSERC) Grant No. ODGP0110050, NATEQ-Quebec Grant No. 02ER3218, and a Canada Research
Chair to the second author, as well as JSPS Grant-In-Aid No. 14654021 to the third author. The
first author benefitted from NSERC and NATEQ scholarships.
Authors’ addresses: F. Panneton and P. L’Ecuyer, Département d’Informatique et de RechercheOpérationnelle, Université de Montréal, C.P. 6128, Succ. Centre-Ville, Montréal, H3C 3J7, Canada;email: {panneton,lecuyer}@iro.umontreal.ca; M. Matsumoto, Department of Mathematics, Graduate School of Science, Hiroshima University, 1-3-1, Kagamiyama, Higashi, Hiroshima, 739-8526,
Japan.
Permission to make digital or hard copies of part or all of this work for personal or classroom use is
granted without fee provided that copies are not made or distributed for profit or direct commercial
advantage and that copies show this notice on the first page or initial screen of a display alongwith the full citation. Copyrights for components of this work owned by others than ACM must be
honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers,
to redistribute to lists, or to use any component of this work in other works requires prior specificpermission and/or a fee. Permissions may be requested from Publications Dept., ACM, Inc., 1515
Broadway, New York, NY 10036 USA, fax: +1 (212) 869-0481, or permissions@acm.org.
C 2006 ACM 0098-3500/06/0300-0001 $5.00
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006, Pages 1–16.

2•

F. Panneton et al.

integer. The main advantages of these generators are that fast implementations
are available and that their mathematical properties can be studied in detailfrom a theoretical perspective. In particular, large-period linear RNGs are easy
to design and the geometrical structure of the set t of all vectors of t successive values produced by the generator, from all its possible initial states, canbe precisely assessed without generating the points explicitly [L’Ecuyer 2004;
Tezuka 1995]. In this article, we are concerned with RNGs based on linear recurrences in F2 , the finite field with two elements, 0 and 1, in which the arithmetic
operations are equivalent to arithmetic modulo 2. Their output sequence is designed to imitate i.i.d. random variables uniformly distributed over the realinterval [0, 1]. Such recurrences are attractive because their implementation
can exploit the binary nature of computers. The set t in this context is a finite
subset of the unit hypercube [0, 1]t and it is customary to require that t covers this unit hypercube very evenly for all t up to some large integer (as large
as possible), because t can be viewed as the “sample space” from which thet-dimensional points are drawn instead of drawing them uniformly from [0, 1]t[L’Ecuyer 2004].
A general framework for representing linear generators over F2 is given by
the matrix linear recurrencexi = Axi−1 ,(1)yi = Bxi ,wui =y i,−1 2− = . y i,0 y i,1 y i,2 · · · ,

(2)
(3)

=1

where xi = (xi,0 , . . . , xi,k−1 )T ∈ Fk2 and yi = ( y i,0 , . . . , y i,w−1 )T ∈ Fw
2 are the k-bit
state and the w-bit output vector at step i, A and B are a k × k transition matrix
and a w × k output transformation matrix (both with elements in F2 ), k and w
are positive integers, and ui ∈ [0, 1) is the output at step i. All operations in (1)and (2) are performed in F2 and each element of F2 is represented as one bit. In
this setup, we havet = (u0 , u1 , . . . , ut−1 ) : x0 ∈ Fk2 .
By making appropriate choices of A and B, several well-known types of generators can be obtained as special cases of this general class, including the
Tausworthe, linear feedback shift register (LFSR), generalized feedback shift
register (GFSR), twisted GFSR, Mersenne twister, and linear cellular automata[L’Ecuyer and Panneton 2002; Panneton 2004].
We briefly recall some (important) well-known facts about this class of RNGs.
The characteristic polynomial of the matrix A can be written as
P (z) = det(A − zI) = z k − α1 z k−1 − · · · − αk−1 z − αk ,where I is the identity matrix and each α j is in F2 . For each j , the sequences{xi, j , i ≥ 0} and { y i, j , i ≥ 0} both obey the linear recurrencexi, j = (α1 xi−1, j + · · · + αk xi−k, j ) mod 2
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

(4)

Generators Based on Linear Recurrences Modulo 2•

3[Niederreiter 1992; L’Ecuyer 1994] (the output sequence { y i, j , i ≥ 0} could alsoobey a recurrence of shorter order in certain cases, depending on B). Here,
we assume that αk = 1, in which case (4) has order k and is purely periodic.
The period length of the recurrence (4) has the upper bound 2k − 1, which isachieved if and only if P (z) is a primitive polynomial over F2 [Niederreiter 1992;
Knuth 1998]. In this article, we are only interested in generators having thismaximal-period property.
The general goal when constructing this type of generator is to find matrices
A and B such that (i) the corresponding recurrence can be implemented very
efficiently by exploiting the special structure of these matrices, (ii) the set t
has excellent uniformity for all t up to some large integer, and (iii) the generator
passes reasonable empirical statistical tests. Typically, the matrix A represents
simple, inexpensive operations such as shifts, rotations, masks, logical ands,
and exclusive-ors (xors) between blocks of bits in the vector xi . The fastest
generators have only a small number of these operations, which implies that
there are very few changes to the bits of the state xi from one step to the next,
and this can often be detected by statistical tests. In particular, if xi contains
many more 0’s than 1’s, the same is likely to hold for xi+1 . To obtain good androbust generators, the matrices A and B in (1) and (2) must perform enough bit
transformations. Therefore, a compromise must be made between speed andquality. The uniformity of t is often measured in terms of its equidistribution
properties, whose definitions are recalled in Section 2.
Our aim in this article is to build a new class of linear generators overF2 that reach a compromise in this sense. The proposed generators perform
more bit transformations and are better equidistributed than (for example) the
Mersenne twister, while having the same period length and approximately thesame speed. To achieve this, we construct a matrix A that implements bit shifts,
bit masks, ands, and xors cleverly spread out in the matrix. For B, in most cases
we simply use Iw×k , the w×w identity matrix to which we append k −w columnsof zeros. We prefer to put all (or most of) the transformations in A because thesetransformations carry over to the following states, whereas those in B do not.
In the next section of the article, we recall basic definitions related to the(standard) use of equidistribution to measure the uniformity of t for linear
generators over F2 . We also specify the uniformity measures adopted in this article. In Section 3, we mention other linear generators over F2 that are alreadyavailable, including combined LFSRs [L’Ecuyer 1999b], GFSRs, the Mersenne
twister [Matsumoto and Nishimura 1998; Nishimura 2000], and combined
twisted GFSRs [L’Ecuyer and Panneton 2002], and we point out the drawbacks
of some of them. Section 4 describes the design of our new generators. Section 5gives specific instances found by a computer search, using REGPOLY [L’Ecuyer
and Panneton 2002], for well-equidistributed generators having state spaces ofvarious sizes. Section 6 deals with implementation issues and speed. In Section 7, we compare the TT800 generator of Matsumoto and Kurita [1994] andthe Mersenne twister of Matsumoto and Nishimura [1998] with two of our new
generators having the same period lengths (2800 −1 and 219937 −1, respectively),in terms of their ability to quickly exit a region of the state space where thefraction of zeros in the state (or in the output vector yi ) is far away from 1/2.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

4•

F. Panneton et al.

The new generators perform much better in these experiments. That is, their
behavior is more in line with randomness and chaos.
2. EQUIDISTRIBUTION AND MEASURES OF QUALITY
A convenient and rather standard way of measuring the uniformity of t for linear RNGs over F2 is as follows. Recall that t has cardinality 2k . If we divide theinterval [0, 1) into 2 equal segments for some positive integer , this determinesa partition of the unit hypercube [0, 1)t into 2t cubic cells of equal size, calleda (t, )-equidissection in base 2. The set t is said to be (t, )-equidistributed if
each cell contains exactly 2k−t of its points. Of course, this is possible only ifdef ≤ ∗t = min(w, k/t).
When the equidistribution holds for all pairs (t, ) that satisfy this condition,
we say that t (and the RNG) is maximally-equidistributed (ME). ME generators have the best possible equidistribution properties in terms of cubic
equidissections. A major motivation for using this type of uniformity measure
for linear generators over F2 is that it can be efficiently computed in this casewithout generating the points explicitly. The idea is simple: The first  bits ofui , . . . , ui+t−1 form a t-bit vector that can be written as Mt, xi for some t × kbinary matrix Mt, because each output bit can be expressed as a linear combination of the bits of the current state, and the (t, )-equidistribution holds if
and only if the matrix Mt, has full rank [Fushimi and Tezuka 1983; L’Ecuyer
1996]. When k is very large this matrix becomes expensive to handle, but inthat case the equidistribution can be verified by a more efficient method (basedon the computation of the shortest nonzero vector in a lattice of formal series),
as explained in Couture and L’Ecuyer [2000]. Large-period ME (or almost ME)
generators have been proposed, for example, by L’Ecuyer [1999b], L’Ecuyer andPanneton [2002], and Panneton and L’Ecuyer [2004]. The aim of this article isto propose new ones with very large periods, which are also faster than thosealready available with comparable period lengths.
For nonME generators, we denote t as the largest dimension t for which tis (t, )-equidistributed, and define the dimension gap for  bits of resolution asδ = t∗ − t ,where t∗ = k/ is an upper bound on the best possible value of t . As measures of uniformity, we consider the worst-case dimension gap and the sum of
dimension gaps, defined as
∞ = max δ
1≤≤w

and

1 =

w


δ .

=1

Aside from equidistribution, it has been strongly advocated that good linear
generators over F2 must have characteristic polynomials P (z) whose numberof nonzero coefficients is not too far from half the degree (i.e., in the vicinity ofk/2) [Compagner 1991; Wang and Compagner 1993]. In particular, generators
for which P (z) is a trinomial or a pentanomial, which have been widely usedin the past, do not satisfy this condition and have been shown to fail rather
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

5

simple statistical tests [Lindholm 1968; Matsumoto and Kurita 1996]. So, asa secondary quality criterion we look at the number of nonzero coefficients in
P (z), which we denote by N1 .
3. OTHER LINEAR GENERATORS OVER F2LFSR generators were studied a long time ago by Tausworthe [1965]. For atrinomial-based LFSR generator, we have w = k, B = I, and A = Ãs (in F2 )for some small positive integer s, where Ã is the matrix that implements therecurrence for which xi,0 = xi−1,q ⊕ xi−1,k−1 and xi, = xi−1,−1 for 0 <  <k, where 0 < q < k and ⊕ represents the “exclusive-or” (xor) operation. A
trinomial-based GFSR generator is defined by the recurrence vi,0 = vi−1,q ⊕
vi−1,r−1 and vi, = vi−1,−1 for 0 <  < r, where k = wr for some integer r, xi
TT
is decomposed into w-bit blocks as xi = (vi,0, . . . , vi,r−1)T , and B = I. Its largestrpossible period length is 2 − 1, which is much smaller than 2k .
A ME trinomial-based LFSR has been proposed already by Tootill et al.
[1973], who introduced the ME notion under the name of asymptotically random. Very fast algorithms are available for implementing trinomial-based
LFSR generators whose parameters satisfy certain conditions, as well as fortrinomial-based and pentanomial-based GFSR generators. However, these generators have far too few nonzero coefficients in their characteristic polynomials(N1 is 3 or 5, which is too small). One way of getting fast generators with alarge value of N1 , proposed by Tezuka and L’Ecuyer [1991] and Wang and Compagner [1993], and pursued by L’Ecuyer [1996, 1999b], is to combine several
trinomial-based LFSR generators of relatively prime period lengths, by bitwise
xor. This gives another LFSR whose characteristic polynomial P (z) is the product of the characteristic polynomials of the components so that it may containmany more nonzero coefficients (e.g., up to 3 J if we combine J trinomial-based
LFSRs and 3 J < k). The polynomial P (z) cannot be primitive in this context,
but the period length can nevertheless be very close to 2k where k is the degreeof P (z). This method is quite effective to build generators with ∞ = 1 = 0
and values of N1 up to a few hundreds on 32-bit computers, but for values inthe thousands or tens of thousands, one needs a large number J of components
and this makes the implementation slower.
The twisted GFSR and Mersenne twister [Matsumoto and Kurita 1994; Matsumoto and Nishimura 1998; Nishimura 2000] provide very efficient implementations of linear generators over F2 with primitive characteristic polynomialsof very large degree k. They are based on a recurrence, similar to that of the
GFSR but slightly more complicated, and a matrix B = I. However, their valuesof N1 are typically much smaller than k/2 and they often have a large value of1 due to the fact that their equidistribution is far from optimal in large dimensions. For example, MT19937 has k = 19937, N1 = 135, and 1 = 6750. Thiscan be improved by combining several twisted GFSRs or Mersenne twisters,
as in L’Ecuyer and Panneton [2002], but at the expense of getting a slower
generator.
Our goal in this article is to build F2 -linear generators with primitive characteristic polynomials, with speed and period length comparable to the Mersenne
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

6•

F. Panneton et al.

Fig. 1. The WELL algorithm.

twister, and for which N1 ≈ k/2 and 1 = 0 (or nearly). We attach the acronym
WELL, for “well-equidistributed long-period linear,” to these new generators.
4. STRUCTURE OF THE NEWLY PROPOSED GENERATORSThe general recurrence for the WELL generators is defined by the algorithm
given in Figure 1, with the following notation. We decompose k as k = rw − pwhere r and p are the unique integers such that r > 0 and 0 ≤ p < w. The
T
T
state vector xi is decomposed into w-bit blocks as xi = (vi,0, . . . , vi,r−1)T where
the last p bits of vi,r−1 are always zero, and T0 , . . . , T7 are w×w binary matricesthat apply linear transformations to these w-bit blocks. The integers k, p, m1 ,m2 , m3 , where 0 < m1 , m2 , m3 < r, and the matrices T0 , . . . , T7 are chosen
so that P (z), the characteristic polynomial of A with degree k = rw − p, isprimitive over F2 . Here, m p is a bit mask that keeps the first w − p bits andsets all p other bits to zero, whereas m̃ p denotes its bitwise complement thatsets the first p bits to zero. The operators “⊕,” “&,” and “←” are the bitwise xor,
the bitwise and, and the assignment statement, respectively. The temporary
variables z0 , . . . , z4 are w-bit vectors.
Let Ti, j,k = Ti Tk ⊕ T j Tk ,




0 0
Iw− p 0Up =and Lw− p =0 0
0 Ip
be matrices of sizes w × w, where I p and Iw− p are the the identity matrices
of size p and w − p, respectively. The matrix A that corresponds to the WELL
algorithm has the following block structure, with six w×w nonzero submatrices
on the first line and four on the second line:⎛⎞
T5,7,0 0 0 . . . T5,7,1 0 . . . T6,7,2 0 . . . T6,7,3 . . . T4 U p T4 Lw− p⎜ T0 0 0 . . . T 1 0 . . . T 2 0 . . . T 3 . . . 0⎟0⎜⎟⎜ 0 I 0 ... 0 0 ... 0 0 ... 0 ... 0⎟0⎜⎟
A = ⎜ 0 0 I ... 0 0 ... 0 0 ... 0 ... 0
⎟.
0
⎜
⎟
⎜
⎟
..
⎝
⎠
.
T
0 0 0 . . . 0 0 . . . 0 0 . . . 0 . . . Lw− p0
As can be deduced from Figure 1, the submatrices T5,7,1 , T6,7,2 , and T6,7,3 arethe (m1 + 1)th, (m2 + 1)th, and (m3 + 1)th w × w blocks of the first line in theblock decomposition of A. If both T0 and T4 have full rank, then A has rank
wr − p.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

7

Table I. Possibilities for the Transformations T j
Transformation Matrix M
M0
M1
M2 (t)−w ≤ t ≤ w
M3 (t)−w ≤ t ≤ w
M4 (a)a ∈ Fw2
M5 (t, b)−w ≤ t ≤ wb ∈ Fw2
M6 (q, t, s, a)
0 ≤ q, s, t < w
a ∈ Fw
2

Implementation of y = Mx = M(x0 , . . . , xk−1 )T
y=0
y=x(x  t) if t ≥ 0y=(x  −t) otherwisex ⊕ (x  t) if t ≥ 0y=x ⊕ (x  −t) otherwise(x  1) ⊕ a if xw−1 = 1y=(x  1)otherwisey=x ⊕ ((x  t)&b) if t ≥ 0
x ⊕ ((x  −t)&b) otherwisey=(((x  q) ⊕ (x  (w − q)))&ds ) ⊕ a if xt = 1(((x  q) ⊕ (x  (w − q)))&ds )otherwise

We take yi = vi,0 as the output vector at step i. Another possibility would
be to take yi = vi,1 , but in our empirical investigations we generally obtained
better equidistribution by taking yi = vi,0 . Taking yi = vi, j for j > 1 would
be equivalent to taking yi = vi,1 , because vi, j = vi−1, j −1 for j = 2, . . . , r − 1,
whereas vi,r is vi−1,r−1 truncated to its first w − p bits.
The transformations (matrices) T j are selected among the possibilities enumerated in Table I. In this table, the bit vector ds has its (s + 1)-th bit
set to zero and all other bits set to one, while “x  t” and “x  t” denote a left shift and a right shift by t bits, respectively. These transformations can be implemented very efficiently by just a few operations on a w-bit
block.
This particular form of A can be motivated intuitively by the following arguments. First, for reasons of efficiency, we want to tap only a small numberof blocks vi−1, j and modify only a small number of blocks vi, j at each iteration.
This is why the majority of rows and columns in the block structure of A are
zero, except for identity matrices in the block subdiagonals whose roles are to
shift the unmodified w-bit blocks by one position. At each step, we use six blocks
to modify the first two, vi,0 and vi,1 . These two blocks will go down the vector andwill be retapped several times in future steps, so their modification has morerepercussions than if we were to modify some other vi, j for j near r. This gives
room for better bit-mixing, which, according to our empirical observations, increases the chances of having a good equidistribution. We use the same blocks
to modify both vi,0 and vi,1 , and even use them in the same way to a certain
extent, again for efficiency reasons. That is, the algorithm first computes z0 ,z1 , and z2 , and then recycles them to obtain z3 and z4 . We tried several othervariants of this general approach with a constraint on the weighted number ofoperations required at each step (giving an appropriate weight to each type ofoperation), and the algorithm retained is the one that gave the best results interms of equidistribution.
The Mersenne twister (MT) generator of Matsumoto and Nishimura [1998]
turns out to be essentially a simplified version of the WELL algorithm, where
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

8•

F. Panneton et al.

T1 = T2 = T3 = T6,7,2 = T6,7,3 = 0, T0 = T5,7,1 = I, and T4 is a companion
matrix whose only nonzero elements are on the first line and on the first subdiagonal (which contains all 1’s). The proposed method is thus a substantial
generalization of MT.
5. SPECIFIC PARAMETERS
Table II lists specific parameters for generators of period length ranging from2512 − 1 to 244497 − 1. The values of the vectors ai are a1 = da442d24, a2 =d3e43ffd, a3 = 8bdcb91e, a4 = 86a9d87e, a5 = a8c296d1, a6 = 5d6b45cc, anda7 = b729fcec. Several of the proposed generators are ME (they have ∞ =
1 = 0). For example, all the generators listed with k = 512, 521, 607, 1024
have this property. We did not find ME generators with B = Iw×k for k =800, 19937, 21701, 23209, 44497, but all those given in Table II have ∞ = 1,
a very small value of 1 (from 1 to 7), and values of N1 not far from k/2. For
comparison, the TT800 generator proposed by Matsumoto and Kurita [1994],with period length 2800 − 1, has 1 = 261 and N1 = 93, whereas the Mersenne
twister MT19937 of Matsumoto and Nishimura [1998] has 1 = 6750 and N1 =135. Table III complements Table II by giving for each generator the values of for which δ > 0.
For the generators that are not ME, it is easy to make them ME by addinga Matsumoto-Kurita tempering to the output. This is done with the following
operations:z ← truncw (xi )z ← z ⊕ ((z  7) & b)yi ← z ⊕ ((z  15) & c)where b and c are carefully selected w-bit vectors and “truncw ” means truncation to the first w bits. The output value ui is generated via equation (3) with thisyi . Adding this tempering to WELL44497a with b = 93dd1400 and c = fa118000
provides the largest ME generator with 32 bits of accuracy known so far. Adding
this tempering to WELL19937a with b = e46e1700 and c = 9b868000 transforms it into the ME generator WELL19937c.
To find the generators listed in Table II, we did a random search using
REGPOLY [L’Ecuyer and Panneton 2002]. It was designed to look for full period generators that use a small number of binary operations. For example, itwould not consider the case where T0 = · · · = T7 = M6 because M6 is the mostexpensive of the matrices of Table I. Once a full-period generator is found, its
equidistribution properties are verified.
To find a full-period generator, we must find a transition matrix having aprimitive characteristic polynomial. To do this, we find the minimal polynomial
over F2 of the sequence { y n,0 }n≥0 generated by the first bit of the yi ’s. If this polynomial is of degree k, then it is the characteristic polynomial of the transition
matrix, otherwise it is a divisor of it. In the former case, we test this polynomial
for primitivity using an algorithm proposed by Rieke et al. [1998] if 2k − 1 is notprime. If 2k −1 is a Mersenne prime, then we can switch to an irreducibility test
which is simpler and equivalent. For this purpose, we used a combination of the
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

9

Table II. Specific Well-Equidistributed Generators
T2
T6

T3
T7

1
N1

k = 512, w = 32, r = 16, p = 0
WELL512aM3 (−16)
M3 (−15)
M3 (−18)1395M3 (−2)M3 (11)M2 (−28)M0
M5 (−5, a1 )0225

k = 521, w = 32, r = 17, p = 23
WELL521aM3 (−13)
M3 (−15)131110
M3 (−13)
M2 (1)
WELL521b
M3 (−21)
M3 (6)
M2 (−10)11107M3 (13)M1
M0
M0
M2 (−5)M2 (−21)
M3 (11)M3 (−13)
M3 (13)

0
265
0
245

k = 607, w = 32, r = 19, p = 1
WELL607a
M3 (19)M3 (11)M1
16
15
14
M3 (18)
WELL607b
M3 (−18)
M3 (−14)
M3 (5)16813
M3 (−24)

M3 (−14)M0
M0
M3 (−1)M1
M3 (−5)M3 (18)
M0

0
295
0
313

k = 800, w = 32, r = 25, p = 0
WELL800aM1
M3 (−15)
M2 (20)141817
M3 (16)
WELL800bM3 (−29)
M2 (−14)
M3 (10)9422M1

M3 (10)M1
M1
M4 (a2 )M3 (−11)
M3 (−28)
M2 (19)M3 (−25)

3
303
3
409

k = 1024, w = 32, r = 32, p = 0
WELL1024a
M1
M3 (8)
M3 (−7)324
10
M3 (−11)
WELL1024b
M3 (−21)
M3 (17)2225
26
M3 (−14)
M3 (−21)

M3 (−19)
M3 (−13)
M4 (a3 )M1

M3 (−14)M0
M3 (15)
M0

0
407
0
475

M3 (1)
M3 (21)M3 (−10)
M3 (−10)

4
8585
5
9679
0
8585

m1

m2

m3

T0
T4

T1
T5

k = 19937, w = 32, r = 624, p = 31
WELL19937a
M3 (−25)
M3 (27)
M2 (9)
70
179 449
M1
M3 (−9)M3 (−21)
WELL19937b
M3 (7)M1
M3 (12)203 613 123
M3 (−19)
M2 (−11) M3 (4)WELL19937c
WELL19937a with M-K temperingb=e46e1700, c=9b868000
k = 21701, w = 32, r = 679, p = 27
WELL21701aM1
M3 (−26)
M3 (−11)151 327 84
M3 (27)M3 (19)M6 (15, 10, 27, a4 )M0
M3 (−16)17609

k = 23209, w = 32, r = 726, p = 23
WELL23209a
M3 (28)M1
M3 (−17)667 43462
M3 (21)
WELL23209b
M4 (a5 )M1
M1
610 175 662
M3 (−26)

M3 (18)M3 (−28)
M6 (15, 30, 15, a6 )
M0

M3 (3)
M3 (−1)M3 (−24)
M3 (16)

3
10871
3
10651

M2 (−26)
M1

7
16883
0
16883

k = 44497, w = 32, r = 1391, p = 15
WELL44497a
M3 (−24)
M3 (30)M3 (−10)23481 229M1
M3 (20)M6 (9, 14, 5, a7 )WELL44497b
WELL44497a with M-K tempering
b=93dd1400, c=fa118000

sieving algorithm described by Brent et al. [2003] and the Berlekamp [1970]
algorithm implemented in the software package ZEN [Chabaud and Lercier
2000]. The method used to compute the equidistribution is the one introduced
by Couture and L’Ecuyer [2000].
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

10•

F. Panneton et al.
Table III. Nonzero Dimension Gaps for the
Proposed WELL’s
WELL
800a
800b
19937a19937b21701a23209a23209b44497a

∞
1
1
1
1
1
1
1
1

1
3
3
4
5
1
3
3
7

{ : δ = 1}
{20, 25, 32}{5, 17, 25}{2, 7, 15, 28}{3, 9, 14, 16, 32}{20}{6, 23, 24}{3, 4, 12}{2, 3, 4, 8, 16, 24, 27}

Fig. 2. An implementation of WELL1024a in C.

6. IMPLEMENTATION AND PERFORMANCE
Figure 2 gives an implementation in C of WELL1024a, whose parameters are
given in Table II. This implementation uses an array of r = 32 unsigned int’s
to store the vi, j vectors. The integer state i is equal to i mod r. Since r is a
power of 2 in this particular case, the operation “modulo 2r ” can be implemented
efficiently using a bit mask, namely, 0x0000001f. The vector vi, j is stored in
STATE[(r − i + j ) mod r]. By decrementing state i at each function call, we
automatically perform the “for” loop of the generator algorithm, saving costly
assignment operations.
Implementations of WELL512a, WELL19937a, WELL19937c, WELL44497a, and WELL44497b are available on the web pages of the first two authors.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

11

Table IV. Time to Generate and
Sum 109 Random Numbers
Generator
WELL512a
WELL1024a
WELL19937a
WELL19937c
WELL44497a
WELL44497b
F2wLFSR2 31 800
F2wLFSR3 7 800
TT800
MT19937
MRG32k3aTime (s)
35.8
35.7
37.1
37.2
40.9
38.8
39.6
35.2
42.7
30.9
97.0

The C code of the last four generators is significantly longer and more complicated than that of WELL512a and WELL1024a, because the order k is nolonger a power of 2, and we are using special tricks to avoid performing several
costly “mod r” operations.
In Table IV, we give the CPU time taken by a program that produces
109 random numbers and sums them using different WELL generatorsfrom Table II: the Mersenne Twister [Matsumoto and Nishimura 1998],
the TT800 [Matsumoto and Kurita 1994], the F2wLFSR3 7 800, and the
F2wLFSR2 31 800 [Panneton and L’Ecuyer 2004], which are all fast generators based on linear recurrences modulo 2. For comparison, we also include
the widely used MRG32k3a [L’Ecuyer 1999a]. For the last five generators, weused the implementations published in the original articles. The test was performed on a computer equipped with a 2.8Ghz Intel Pentium 4 processor and a
Linux operating system, and the program was compiled using gcc with the -O2
optimization flag. We repeated the timing experiments five times and the results agreed to within 0.5 seconds of accuracy. The numbers are rounded to thenearest integer. Interestingly, the WELL44497b generator (with tempering) isslightly faster than its basic version WELL44497a (without tempering). This
strange behavior is probably an artefact of the compiler optimization methods, which operate differently in the two cases. The biggest WELL generators
are almost as fast as the smaller ones. On the other hand, they use a larger
amount of memory. So, it is not necessarily true that WELL44497b should always be preferred over WELL1024, for example, even if it has a much longer
period.
The WELL generators mentioned in Table IV successfully passed all thestatistical tests included in the batteries Smallcrush, Crush, and Bigcrush of
TestU01 [L’Ecuyer and Simard 2002], except those that look for linear dependencies in a long sequence of bits, such as the matrix-rank test [Marsaglia
1985] for very large binary matrices and the linear complexity tests [Erdmann
1992]. This is, in fact, a limitation of all F2 -linear generators, including theMersenne twister, the TT800, etc. Because of their linear nature, the sequences
produced by these generators just cannot have the linear complexity of a truly
random sequence. This is definitely unacceptable in cryptology, for example,
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

12•

F. Panneton et al.

Fig. 3. γn,100 for the WELL800a (top line) and TT800 (lower line).

but is quite acceptable for the vast majority of simulation applications if thelinear dependencies are of long range and high order.
7. ESCAPING ZEROLAND
In this section, we examine how quickly the WELL generators can escape froma bad initialization (e.g., one that contains only a few bits set to 1). We comparethem with other generators for which the initialization often has a long lasting
impact, in the sense that the fraction of bits set to 1 in both the state and theoutput remains small for a large number of steps if the initial number is verysmall.
Our experiments are based on the following framework. Let S1 = {e j : 1 ≤j ≤ k}, where e j denotes the j -th unit vector. This is the set of states for which( j)a single bit is set to 1. Let yi be the output vector yi at step i when the initial( j)state is x0 = e j , and let ui be the corresponding output value ui . We denoteby H(x) the Hamming weight of a bit vector x (i.e., the number of its bits thatare set to 1) and define
γn, p =

n+
p−1 
k

 ( j )1
H yi .
pkw i=n j =1

(5)This moving average represents the fraction of bits that are 1 in the p successive
output vectors yn , . . . , yn+ p−1 , averaged over all initial seeds in S1 . Under thenull hypothesis that the yi ’s are independent bit vectors uniformly distributed
over Fw
2 , γn, p should be approximately normally distributed with mean 1/2 andvariance 1/(4 pkw).
Figure 3 illustrates the behavior of γn,100 for n = 100, . . . , 105 , for the
WELL800a and TT800 generators. We see that WELL800a reaches the neighborhood of 0.5 faster than TT800. For the latter, γn,100 is often below 0.49 for n
up to more than 60000. Here, the standard deviation of γn,100 should be 1/3200,
so 0.49 is 32 standard deviations away from the mean.
In Figure 4, we see a similar comparison between WELL19937a andMT19937, with γn,1000 and n up to one million. Again, the WELL generator
is much better behaved. With MT19937, γn,1000 has an easily visible negative
bias for up to, at least, n = 700000 steps.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

13

Fig. 4. γn,1000 for the WELL19937a (top line) and MT19937 (lower line).

Fig. 5. γn,5 for the WELL800a (top line) and WELL19937a (lower line).

Although not visible in the scale of these two figures, the WELL generators
also need a certain number of steps to bring their states close to the “center” ofthe state space, where roughly half of the bits are 1’s. To show this, Figure 5 gives
a close-up view of the early behavior of γn,5 for WELL800a and WELL19937a.
There is still a visible (and significant) amount of dependence from the initialstate, but the number of steps required to clear up the effect is approximately
1000 times smaller than for TT800 and MT19937.
( j)
( j)
We also computed similar figures, but with ui instead of H(yi )/w in (5)(i.e., the average output value instead of the overall fraction of bits set to 1),
and the results were almost identical. That alternative gives more weight tothe more significant bits.
The behavior of the Hamming weight that we just examined is closely related
to the diffusion capacity of generators. Crudely speaking, a generator is said tohave good diffusion capacity if it produces two sequences that are very different
from each other when started from two seeds that are very close to each other[Shannon 1949]. A precise definition depends on how we measure the closenessof seeds and the differences of sequences.
For the generators considered in this article, we may define the distance
between two seeds or two output vectors as the number of bits that differ between them. This is the Hamming distance. If two seeds x0 and x0
are at distance 1, then x0 ⊕ x0 = e j for some j . If yi and yi are the corresponding output vectors at step i, then because of the linearity, we have
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

14•

F. Panneton et al.

yi ⊕ yi = BAi (x0 ⊕ x0 ) = BAi e j = yi , so the distance between yi and yi is the( j)Hamming weight of yi . In this sense, the average Hamming weight definedin (5) also measures the diffusion capacity of the generator.
What is the potential impact of this weakness of MT and TGFSR generators on simulation results for practical applications? If the generator’s state is
initialized with k independent, uniform (truly) random bits, the probability ofhitting a bad state is probably quite small (a more precise statement would require a formal definition of a “bad state” and further analysis). But chances still
exist that the generator eventually hits a bad state during a simulation. Moreover, such huge generators are rarely seeded by k truly random bits in practice,and seeds that contain many zeros are perhaps more likely to be used because
they are easier to write down. This could lead to totally wrong simulation results. The WELL generators provide a substantial improvement in this aspect.
( j)

8. CONCLUSIONS
Random number generators based on linear recurrences modulo 2 are popularand very convenient for simulation, mainly because of their high speed andbecause the quality of specific instances can be studied in detail from the mathematical viewpoint. Other generators from this family have been proposed in
the past, for instance, the Mersenne twister and the combined LFSRs. Some ofthem have huge period lengths, so their states must be represented over a large
number of bits, and require only a few operations on 32-bit (or 64-bit) words
to go from one state to the next. This means that only a small fraction of thebits in the state are modified at each step. In this article, we have proposed aclass of generators that do better than the previous ones in this aspect, withoutcompromising the period length and speed. The improvement was assessed bya measure of equidistribution, by counting the fraction of nonzero coefficientsin the characteristic polynomial of the recurrence, and by an empirical comparison of the diffusion capacity of the old and new generators. Of course, like all
linear generators over F2 , these generators cannot pass a statistical test thatmeasures the linear complexity of a long sequence of bits that they produce. But
they are, nevertheless, very useful and safe for the vast majority of simulation
applications, especially when the generator’s speed is important.
ACKNOWLEDGMENTS

We thank Takuji Nishimura, who gave us useful ideas to improve the performance of our algorithms for checking primitivity and equidistribution, Étienne
Marcotte, who produced the figures, Hiroshi Haramoto, who found some mistakes in an earlier version, and the anonymous reviewers whose comments
helped improve the article.
REFERENCES
BERLEKAMP, E. R. 1970. Factoring polynomials over large finite fields. Math. Comput. 24, 713–735.
BRENT, R. P., LARVALA, S., AND ZIMMERMANN, P. 2003. A fast algorithm for testing reducibility oftrinomials mod 2 and some new primitive trinomials of degree 3021377. Math. Comput. 72, 243,
1443–1452.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

Generators Based on Linear Recurrences Modulo 2•

15

CHABAUD, F. AND LERCIER, R. 2000. A toolbox for fast computation in finite extension over finite
rings. Software User’s Guide. See http://zenfact.sourceforge.net/.
COMPAGNER, A. 1991. The hierarchy of correlations in random binary sequences. J. Stat. Phy. 63,
883–896.
COUTURE, R. AND L’ECUYER, P. 2000. Lattice computations for random numbers. Math. Comput. 69, 230, 757–765.
ERDMANN, E. D. 1992. Empirical tests of binary keystreams. M.S. thesis, Department of Mathematics, Royal Holloway and Bedford New College, University of London.
FUSHIMI, M. AND TEZUKA, S. 1983. The k-distribution of generalized feedback shift register pseudorandom numbers. Commun. ACM 26, 7, 516–523.
KNUTH, D. E. 1998. The Art of Computer Programming, vol. 2: Seminumerical Algorithms, 3rd
ed. Addison-Wesley, Reading, Mass.
L’ECUYER, P. 1994. Uniform random number generation. Ann. Oper. Res. 53, 77–120.
L’ECUYER, P. 1996. Maximally equidistributed combined Tausworthe generators. Math. Comput. 65, 213, 203–213.
L’ECUYER, P. 1999a. Good parameters and implementations for combined multiple recursive random number generators. Oper. Res. 47, 1, 159–164.
L’ECUYER, P. 1999b. Tables of maximally equidistributed combined LFSR generators. Math. Comput. 68, 225, 261–269.
L’ECUYER, P. 2004. Random number generation. In Handbook of Computational Statistics, J. E.
Gentle et al., eds. Springer Verlag, Berlin, 35–70.
L’ECUYER, P. AND PANNETON, F. 2002. Construction of equidistributed generators based on linear
recurrences modulo 2. In Monte Carlo and Quasi-Monte Carlo Methods 2000, K.-T. Fang et al.,eds. Springer Verlag, Berlin, 318–330.
L’ECUYER, P. AND SIMARD, R. 2002. TestU01: A Software Library in ANSI C for Empirical Testing ofRandom Number Generators. Software User’s Guide. http://www.iro.umontreal.ca/∼lecuyer.
LINDHOLM, J. H. 1968. An analysis of the pseudo-randomness properties of subsequences of long
m-sequences. IEEE Trans. Inf. Theor. IT-14, 4, 569–576.
MARSAGLIA, G. 1985. A current view of random number generators. In Computer Science andStatistics, Sixteenth Symposium on the Interface. Elsevier Science Publishers, North-Holland,
Amsterdam, 3–10.
MATSUMOTO, M. AND KURITA, Y. 1994. Twisted GFSR generators II. ACM Trans. Model. Comput.
Simul. 4, 3, 254–266.
MATSUMOTO, M. AND KURITA, Y. 1996. Strong deviations from randomness in m-sequences basedon trinomials. ACM Trans. Model. Comput. Simul. 6, 2, 99–106.
MATSUMOTO, M. AND NISHIMURA, T. 1998. Mersenne twister: A 623-dimensionally equidistributed
uniform pseudo-random number generator. ACM Trans. Model. Comput. Simul. 8, 1, 3–
30.
NIEDERREITER, H. 1992. Random Number Generation and Quasi-Monte Carlo Methods. SIAM
CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 63. SIAM, Philadelphia.
NISHIMURA, T. 2000. Tables of 64-bit Mersenne twisters. ACM Trans. Model. Comput. Simul. 10, 4,
348–357.
PANNETON, F. 2004. Construction d’ensembles de points basée sur des récurrences linéaires dans
un corps fini de caractéristique 2 pour la simulation Monte Carlo et l’intégration quasi-Monte
Carlo. Ph.D. thesis, Département d’informatique et de recherche opérationnelle, Université deMontréal, Canada.
PANNETON, F. AND L’ECUYER, P. 2004. Random number generators based on linear recurrences in
F2w . In Monte Carlo and Quasi-Monte Carlo Methods 2002, H. Niederreiter, ed. Springer Verlag,
Berlin, 367–378.
RIEKE, A., SADEGHI, A.-R., AND POGUNTKE, W. 1998. On primitivity tests for polynomials. In Proceedings of the 1998 IEEE International Symposium on Information Theory. Cambridge, Mass.
SHANNON, C. E. 1949. Communication theory of secrecy systems. Bell Syst. Tech. J. 28, 656–715.
TAUSWORTHE, R. C. 1965. Random numbers generated by linear recurrence modulo two. Math.
Comput. 19, 201–209.
TEZUKA, S. 1995. Uniform Random Numbers: Theory and Practice. Kluwer Academic, Norwell,
Mass.
ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.

16•

F. Panneton et al.

TEZUKA, S. AND L’ECUYER, P. 1991. Efficient and portable combined Tausworthe random number
generators. ACM Trans. Model. Comput. Simul. 1, 2, 99–112.
TOOTILL, J. P. R., ROBINSON, W. D., AND EAGLE, D. J. 1973. An asymptotically random Tausworthesequence. J. ACM 20, 469–481.
WANG, D. AND COMPAGNER, A. 1993. On the use of reducible polynomials as random number generators. Math. Comput. 60, 363–374.
Received July 2004; revised April 2005; accepted April 2005

ACM Transactions on Mathematical Software, Vol. 32, No. 1, March 2006.