On the Xorshift Random Number Generators
FRANÇOIS PANNETON and PIERRE L’ECUYER
Université de Montréal

G. Marsaglia introduced recently a class of very fast xorshift random number generators, whose
implementation uses three “xorshift” operations. They belong to a large family of generators basedon linear recurrences modulo 2, which also includes shift-register generators, the Mersenne twister,and several others. In this paper, we analyze the theoretical properties of xorshift generators,
search for the best ones with respect to the equidistribution criterion, and test them empirically.
We find that the vast majority of xorshift generators with only three xorshift operations, including
those having good equidistribution, fail several simple statistical tests. We also discuss generatorswith more than three xorshifts.
Categories and Subject Descriptors: G.4 [Mathematical Software]: Algorithm design andanalysis; I.6 [Computing Methodologies]: Simulation and Modeling
General Terms: Algorithms
Additional Key Words and Phrases: Random number generation, xorshift, linear feedback shift
register, linear recurrence modulo 2

1.

INTRODUCTION

Marsaglia [2003] proposed a class of very fast uniform random number generators(RNGs) called xorshift. The state of a xorshift generator is a vector of bits. At each
step, the next state is obtained by applying a certain number of xorshift operations
to w-bit blocks in the current state, where w = 32 or 64, and a xorshift operation
is defined as follows: replace the w-bit block by a bitwise xor (exclusive or) of theoriginal block with a shifted copy of itself by a positions either to the right or to
the left, where 0 < a < w.
Xorshifts are linear operations. The left shift of a w-bit vector x by one bit,x  1, can also be written as Lx where L is the w × w matrix with ones on its
main superdiagonal and zeros elsewhere. Similarly, the right shift x  1 can be
written as Rx where R has ones on its main subdiagonal and zeros elsewhere.
Matrices of the forms (I + La ) and (I + Ra ), where a ∈ {1, . . . , w − 1}, are called
left and right xorshift matrices, respectively. They represent left and right a-bit
xorshift operations.
Département d’Informatique et de Recherche Opérationnelle, Université de Montréal, C.P.
6128, Succ. Centre-Ville, Montréal, H3C 3J7, Canada, e-mail: panneton@iro.umontreal.ca,
lecuyer@iro.umontreal.ca
Permission to make digital/hard copy of all or part of this material without fee for personalor classroom use provided that the copies are not made or distributed for profit or commercial
advantage, the ACM copyright/server notice, the title of the publication, and its date appear, andnotice is given that copying is by permission of the ACM, Inc. To copy otherwise, to republish,
to post on servers, or to redistribute to lists requires prior specific permission and/or a fee.
c 20YY ACM 0000-0000/20YY/0000-0001 $5.00
ACM Journal Name, Vol. V, No. N, Month 20YY, Pages 1–16.

2·

F. Panneton and P. L’Ecuyer

A xorshift generator is defined by a recurrence of the form
vi =

p
X

Ãj vi−mj mod 2

(1)

j=1

where p is a positive integer, the vi ’s are w-bit vectors, the mj ’s are positive integers,
and Ãj is either the identity or the product of νj xorshift matrices for some νj ≥ 0,for each j (Ãj is the zero matrix if νj = 0). If we define r = max1≤j≤p mj ,T
the generator’sP
state at step i can be written as xi = (vi−r+1, . . . , viT )T and thew−`T
output is ui = `=1 vi,`−1 2 where vi = (vi,0 , . . . , vi,w−1 ) . The output sequence{ui , i ≥ 0} is supposed to imitate i.i.d. random variables uniformly distributed overthe interval [0, 1]. We can rewrite (1) as
vi =

r
X

Aj vi−j mod 2

j=1

where Aj =

X

Ãl .

(2)

{l:ml =j}This is a special case of the multiple recursive matrix method defined in Niederreiter[1995].
Marsaglia [2003] considers three types of xorshift generators, mostly with ν1 +· · · + νr = 3 (i.e., exactly three xorshift operations). The type I generators haver = 1 and A1 is the product of three xorshift matrices. For example, one may
have A1 = (I + L13 )(I + R17 )(I + L5 ), which Marsaglia says is “one of my favorite
choices.” For the type II generators, r > 1, Ar 6= 0 and only one other matrix
Aj is nonzero. For example, one could have r = 4, A4 = (I + R7 )(I + L11 ), and
A1 = (I + L20 ). The type III generators have r > 1, Ar 6= 0, exactly three matrices
Aj are xorshift matrices, and the others are zero. For example, one may haver = 12, A2 = (I + L7 ), A3 = (I + R11 ), and A12 = (I + L21 ).
These generators are extremely fast and it is easy to find parameter values forwhich they have full period length 2rw − 1. This period length can be made verylarge by selecting a large r and (say) w = 32. But as we all know, a long period
does not suffice to have a high-quality generator. In his paper, Marsaglia says:“Although I have only tested a few of them, any one of the 648 choices above islikely to provide a very fast, simple, high quality RNG.” So, how robust and reliable
are they?
To answer this question, in this paper, we study the xorshift generators fromboth the theoretical and empirical perspectives. These generators form a subclassof a large and well-known family of RNGs based on linear recurrences modulo 2,
as described in Section 2. In Section 3, we recall standard equidistribution criteriathat are widely used to measure the uniformity and independence of these linear
generators. In Section 4, we establish some properties of xorshift matrices andgenerators. In particular, we identify classes of generators having the same period
length, others having the same equidistribution, and others that are completely
equivalent. In Section 5, we assess the equidistribution of specific instances proposed by Marsaglia and submit them to statistical tests. Many of them have very
bad equidistribution and fail some tests spectacularly. In Section 6, we perform asearch for the best generators of type I, II, and III with respect to equidistribution
criteria, under the constraint that they use only three xorshift operations. We test
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·3

the best ones and find that they are statistically weak. In Section 7, we briefly lookat xorshift generators with more than three xorshifts.
2.

MATRIX LINEAR RECURRENCES MODULO 2

The xorshift generators belong to a large class of RNGs based on the following typeof linear recurrences modulo 2:
xi = Axi−1 mod 2,
yi = Bxi mod 2,wXui =yi,`−1 2−` = .yi,0 yi,1 yi,2 · · · ,

(3)
(4)
(5)

`=1

where xi = (xi,0 , . . . , xi,k−1 )T and yi = (yi,0 , . . . , yi,w−1 )T are the k-bit state andthe w-bit output vector at step i, A is a k × k binary transition matrix, B is a w × kbinary output transformation matrix, k and w are positive integers, and ui ∈ [0, 1) is
the output at step i. By appropriate choices of A and B, several types of generators
can be obtained as special cases of this general class; for instance the Tausworthe,linear feedback shift register (LFSR), generalized feedback shift register (GFSR),
twisted GFSR, Mersenne twister, WELL, etc. (see, e.g., L’Ecuyer 2004; Panneton
et al. 2005, and Panneton 2004 for the details). Important advantages of thesegenerators are that fast implementations are available (they exploit the binary
nature of computers) and that their mathematical properties are well understoodfrom a theoretical perspective.
Let P (z) = det(A − Iz) be the characteristic polynomial of A. It is well-known
that the recurrence (3) has period length 2k − 1 (its maximal possible value) if andonly if P (z) is a primitive polynomial modulo 2 [Niederreiter 1992; Knuth 1998].
Xorshift generators fit this setting by taking k = rw,


0
I
··· 0
.. 
..
 ...
.
. ,
A=
 0
0
··· I Ar Ar−1 · · · A1
T
xi = (vi−r+1, . . . , viT )T , and yi = vi . The matrix B contains a w × w identity I in
its upper left cornerzeros elsewhere. The characteristic polynomial of this A is
Pandr
P (z) = det(z r I + j=1 z r−j Aj ). While completing this paper, we found a related
article by Brent [2004a], who also reports the similitudes between the xorshift and
LFSR generators.
Note that replacing the identity I in the upper left corner of B by an arbitrary
invertible w × w matrix B̃, i.e., defining yi = B̃vi , is equivalent to replacing each
Aj in (2) by B̃Aj B̃−1 , which gives the recurrence

yi =

r
X

B̃Aj B̃−1 yi−j mod 2.

(6)

j=1

This does not change the characteristic polynomial P (z).
ACM Journal Name, Vol. V, No. N, Month 20YY.

4

·

3.

EQUIDISTRIBUTION

F. Panneton and P. L’Ecuyer

For an arbitrary integer t > 0, a vector of t successive output values of an RNG
would ideally behave as a random vector uniformly distributed over the t-dimensional
unit hypercube [0, 1]t . But for an RNG defined by (3)–(5), all these vectors must
belong to the finite set
Ψt = {(u0 , u1 , . . . , ut−1 ) : x0 ∈ {0, 1}k }.
It is customary to require that Ψt covers [0, 1]t very evenly for all t up to some large
integer, so that drawing a point uniformly from the “sample space” Ψt gives a good
approximation of a uniform random variable over [0, 1]t [L’Ecuyer 2004]. For thisclass of RNGs, a standard way of assessing the uniformity of Ψt is the following.
Recall that Ψt has cardinality 2k . If we divide the interval [0, 1) into 2` equal
segments for some positive integer `, this determines a partition of the unit hypercube [0, 1)t into 2t` cubic cells of equal size, called a (t, `)-equidissection in base 2.
The set Ψt is said to be (t, `)-equidistributed if each cell contains exactly 2k−t` ofits points. This can be verified by expressing the t` bits that determine the cellnumber of a point as a linear combination of the k bits of the state, and checking if the corresponding matrix has full rank [Fushimi and Tezuka 1983; L’Ecuyer
1996]. For a fixed resolution `, let t` denote the largest value of t such that Ψt is
def(t, `)-equidistributed. A theoretical upper bound on t` is t∗` = bk/`c. We define
the dimension gap in resolution ` as δ` = t∗` − 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
X

δ` .

`=1

A small value of ∆1 or ∆∞ indicates good uniformity. A generator is called maximally equidistributed (ME) if ∆1 = 0 [L’Ecuyer 1996]. ME generators have thebest possible equidistribution properties in terms of cubic equidissections.
4.

PROPERTIES OF XORSHIFT GENERATORSIn this section, we examine some theoretical properties of xorshift matrices and show
how these properties affect the period length and the equidistribution of xorshift
generators.
4.1

General Properties

Our first proposition implies that to reach the maximal period 2rw − 1, a xorshift
generator must contain both left and right xorshift matrices.
Proposition 4.1. If the nonzero matrices Ai are products and/or sums of either
all left or all right xorshift matrices, then P (z) cannot be irreducible, and therefore
cannot be primitive.
Proof. If the nonzero Ai ’s are all P
products and/or sums of xorshift matrices onr
the same side, then the matrix z r I + j=1 z r−j Aj is triangular, hence P (z) is the
product of its main diagonal, i.e. a product of z’s and (1 − z)’s.
Consider a generator defined by (3)—(5) where B has an arbitrary w×w matrix B̃
in its upper left corner and zeros elsewhere. Let Ψt (A, B̃) denote the corresponding
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·5

set Ψt . Our next result says that adding a right xorshift to the output of such agenerator preserves its equidistribution, and therefore does not change its values of∆1 and ∆∞ .
Proposition 4.2. For ` ≤ w, the set Ψt (A, (I + Ra )B̃) is (t, `)-equidistributedfor 0 < a < w if and only if Ψt (A, B̃) is.
Proof. The t` bits that determine in which box of the (t, `)-equidissection the
point (u0 , . . . , ut−1 ) will fall are the ` most significant bits of each of its coordinates.
TLet y(t, `) = (y0,0 , . . . , y0,`−1 , . . . , yt−1,0 , . . . , yt−1,`−1 ) be the vector of those bits
and observe that y(t, `) = M(t, `)x0 where M(t, `) is the t`×k binary matrix whoserow `i + j is row j of BAi , for 1 ≤ j ≤ ` and 0 ≤ i ≤ t − 1. The set Ψt is (t, `)equidistributed if and only if each possibility for y(t, `) occurs the same number oftimes when x0 runs through all of its 2k possibilities. This happens if and only if
the matrix M(t, `) has full rank, t`.
If we consider Ψt (A, T̃B̃) instead of Ψt (A, B̃), where T̃ = I + Ra , we must
redefine row `i + j of M(t, `) as row j of TBAi , where T is a matrix with T̃ in
its upper-left corner and zeros elsewhere. This amounts to left-multiplying M(t, `)by a block-diagonal t` × t` matrix D comprised of t diagonal blocks which are all
equal to the ` × ` upper-left corner of T̃. But for any a > 0, the matrix T̃ islower-triangular with ones all over its main diagonal, so D also has this property
and is thus invertible. Therefore, left-multiplying M(t, `) by D does not change its
rank and this completes the proof.
Two matrices M1 and M2 are similar if there exists an invertible matrix C suchthat M1 = CM2 C−1 . Similarity is an equivalence relationship, denoted M1 ∼ M2 .
Similar transition matrices A have the same characteristic polynomial and define
generators with the same period length. The next proposition exploits this fact toderive several useful properties of xorshift matrices.
Proposition 4.3. Let H, H1 , H2 , H3 , K1 , K2 , K3 ∈ {L, R} and such that Hi 6=
Ki for i = 1, 2, 3. We have that:(a) (I + La ) ∼ (I + Ra );(b) (I + Ha )(I + Hb ) = (I + Hb )(I + Ha );(c) (I + Hc3 )(I + Hb2 )(I + Ha1 ) ∼ (I + Ha1 )(I + Hc3 )(I + Hb2 );(d) (I + Hc3 )(I + Hb2 )(I + Ha1 ) ∼ (I + Kc3 )(I + Kb2 )(I + Ka1 );(e) if H1 , H2 and H3 are not all the same matrix,
then (I + Hc3 )(I + Hb2 )(I + Ha1 ) ∼ (I + Hc3 )(I + Ha1 )(I + Hb2 ).
Proof. Let P denote the square matrix with ones on its main antidiagonal andzeros elsewhere, i.e., the identity matrix with its columns in reverse order. We have
(I + La ) = P(I + Ra )P−1 and this implies (a) and (d). Part (b) can be verified bymultiplying the terms on each side and using the fact that Ha Hb = Ha+b = Hb Ha .
For part (c), we have (I + Hc3 )(I + Hb2 )(I + Ha1 ) ∼ (I + Ha1 )(I + Hc3 )(I + Hb2 )(I +Ha1 )(I + Ha1 )−1 ∼ (I + Ha1 )(I + Hc3 )(I + Hb2 ). Part (e) is left as an exercise to thereader.
By using the properties listed in Proposition 4.3, one can easily show the following
proposition.
ACM Journal Name, Vol. V, No. N, Month 20YY.

·6

F. Panneton and P. L’Ecuyer

Proposition 4.4. The eight matrices
X1 = (I + Lc )(I + Rb )(I + La ), X2 = (I + La )(I + Rb )(I + Lc ),X3 = (I + Rc )(I + Lb )(I + Ra ), X4 = (I + Ra )(I + Lb )(I + Rc ),X5 = (I + Rb )(I + Lc )(I + La ), X6 = (I + Lb )(I + Ra )(I + Rc ),X7 = (I + Lc )(I + La )(I + Rb ), X8 = (I + Ra )(I + Rc )(I + Lb ),where a, b, c are in {1, . . . , w − 1}, are all similar.
Proof. Let P be as in the proof of Proposition 4.3 and Cj , 1 ≤ j ≤ 8, be
the matrices such that X1 = Cj Xj C−1
j . Using Proposition 4.3, we show that:
C2 = (I + La )−1 (I + Lc ), C3 = P, C4 = P(I + Ra )−1 (I + Rc ), C5 = (I + Lc ),C6 = (I + Lc )P, C7 = (I + La )−1 and C8 = (I + La )−1 P.
We will make use of this last proposition in our examination of the three typesof xorshift generators proposed by Marsaglia [2003].
4.2

Generators of type I

Because the matrices X1 , . . . , X8 are all similar, by finding a triple (a, b, c) for which
the generator based on X1 has full period, we in fact identify eight full-period
generators. Marsaglia [2003] recognizes this fact (without providing a proof), but
only for the matrices X1 , . . . , X6 , and he includes the matrices (I + Rb )(I + La )(I +Lc ) and (I + Lb )(I + Rc )(I + Ra ), which are in fact a repetition X5 and X6 .
The next proposition implies that to verify the equidistribution of all type I
generators defined by a triple (a, b, c), we only need to verify the equidistributionof the generators based on X1 , X2 , X3 , and X5 .
Proposition 4.5. For a given triple (a, b, c), the type I generators with A1 =X7 and A1 = X5 have exactly the same equidistribution properties, and the type I
generators with A1 equal to either X3 , X4 , X6 , or X8 also have exactly the sameequidistribution properties.
Proof. By (6), taking A1 = X7 and yi = (I + Rb )vi is equivalent to using
the recurrence yi = (I + Rb )X7 (I + Rb )−1 yi−1 = (I + Rb )(I + Lc )(I + La )(I +Rb )(I + Rb )−1 yi−1 = X5 yi−1 . This means that Ψt (X5 , I) = Ψt (X7 , I + Rb ) for
all t. It then follows from Proposition 4.2 that the generators with X7 and X5
have identical equidistribution properties. A similar argument applies to the typeI generators produced by X3 , X4 , X6 , and X8 .
4.3

Generators of type II

A generator of type II is implemented via a special version of (2) which can be
written as:vi = Gvi−m ⊕ Hvi−r(7)where G and H are w×w matrices and ⊕ denotes the bitwise exclusive-or operation.
We denote this generator by the triple (G, H, m). Property (6) implies that for agiven non-singular w × w matrix C, (CGC−1 , CHC−1 , m) provides a full-period
generator if and only if (G, H, m) provides one.
Table I lists all possibilities of G and H. They are grouped in a way thatwithin each group, they all have the same characteristic polynomial through the
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·

7

similarity matrix C. The last column indicates the generators having the sameequidistribution. Marsaglia [2003] only mentions generators for which G = G1 =(I + Rb )(I + La ), H = H1 = (I + Rc ), and m = 1. The following theorem tells usthat it is worthless to search for full-period generators of type II with H = I.

Table I. List of possible combinations (G, H), regrouped in classes of generators having the same
period length.
i
1
2
3
4
5
6
7
8
9
10
11
12
13–14
15–16
17–20

Gi
(I + Rb )(I + La )(I + Lb )(I + Ra )(I + La )(I + Rb )(I + Ra )(I + Lb )(I + Rc )(I + Lc )(I + Rc )(I + Lc )(I + Ra )(I + Rb )(I + La )(I + Lb )(I + Lc )(I + Rc )I
IIHi(I + Rc )(I + Lc )(I + Rc )(I + Lc )(I + Rb )(I + La )(I + Lb )(I + Ra )(I + La )(I + Rb )(I + Ra )(I + Lb )(I + Lc )(I + Rc )(I + Ra )(I + Rb )(I + La )(I + Lb )Xj , j ∈ {1, 2}
Xj , j ∈ {5, 7}
Xj , j ∈ {3, 4, 6, 8}

CIP(I + Rb )(I + Lb )PI
P(I + Rb )(I + Lb )PI
PI
P

Equidistribution
same as i = 3

See Section 4same for i = 15, 16
same for 17 ≤ i ≤ 20

same as i = 1
same as i = 7
same as i = 5

Theorem 4.6. Let G be any binary w × w matrix where w > 1. The recurrence(7) with H = I cannot have period length 2rw − 1.
Proof. For n equal to any power of 2, we denote by Fn the finite field with nelements, Fw
n its w-fold cartesian product, and Fn [z] the space of polynomials withcoefficients in Fn . Let Q(z) = det(G − Iz) be the characteristic polynomial of G
over F2 and suppose that Q(z) is irreducible over F2 . Let η ∈ F2w be a root of
Q(z). The elements 1, η, . . . , η w−1 form a basis of F2w over F2 .
We can show (see the proof of Theorem 2 in Matsumoto and Kurita 1992) thatfor a nonzero vector in t ∈ Fw
2 , the homomorphism
φ : F2 w → Fw2
η l 7→ Gl t,
where 0 6= t ∈ Fw
2 is fixed, is also an isomorphism. By applying the inverse of φ torecurrence (7), we obtain the linear recurrenceφ−1 (xi ) = φ−1 (Gxi−m ) + φ−1 (xi−r ),which can be rewritten as the following linear recurrence in F2w :xi = ηxi−m + xi−r ,(8)where xi = φ−1 (xi ) ∈ F2w , because φ−1 (Gxi−m ) = ηxi−m (see Theorem 2 in Matsumoto and Kurita 1992). Notice that (7) and (8) have the same period length.
ACM Journal Name, Vol. V, No. N, Month 20YY.

8·

F. Panneton and P. L’Ecuyer

The characteristic polynomial over F2w of (8) is
P (z) = z r − ηz r−m − 1 ∈ F2w [z]By Theorem 3.18 of Lidl and Niederreiter [1994], because (−1)r P (0) = (−1)r+1 is
not a primitive element of F2w , the polynomial P (z) is not primitive over F2w . Thisimplies that the period of (8) cannot be 2rw − 1, and similarly for (7).
Now, suppose that Q(z) is not irreducible. Let q = r − m. We know that
P (z) = det(z r I − z q G − I) = det(z q (z m I − G − z −q I))= det(z q I) det((z m − z −q )I − G)= z wq Q(z m − z −q ) ∈ F2 [z]Notice that Q(z m − z −q ) is not in F2 [z], so we do not know yet if P (z) is irreducible
over F2 or not. On the other hand, Q(z m − z −q ) ∈ L2 where L2 is the field offormal Laurent series with coefficients in F2 . Because Q(z) is not irreducible, we
can decompose it in c > 1 factors Qi (z) ∈ F2 [z], each of degree di > 0.
Let h(z) = hn z n + hn−1 z n−1 + · · · ∈ L2 . We define the function p(h(z)) =min{i : hi 6= 0}. For h(z) ∈ L2 to be in F2 [z], it is necessary that p(h(z)) ≥ 0.
Observe that, because Q(z) is of degree w, we have that w = d1 + · · · + dc andp(Qi (z m + z −q )) = −qdi . Let Q̄i (z) = z qdi Qi (z m + z −q ). Then, p(Q̄i (z)) = 0,
which implies that Q̄i (z) ∈ F2 [z]. We can develop
z wq Q(z m − z −q ) = z wqcY

Qi (z m − z −q )

i=1

=
=

c
Y
i=1
c
Y

z qdi Qi (z m − z −q )Q̄i (z).

i=1

The last equality implies that P (z) is not irreducible over F2 because it can be
decomposed in c > 1 factors Q̄i (z) ∈ F2 [z], i = 1, . . . , c. This shows that if Q(z)is not irreducible over F2 , then neither is P (z), and (8) cannot have the maximal
period 2rw − 1 in that case.
So far, we have discussed the choice of H and G, but not the choice of m. For
the case where G = I, conditions on m for the recurrence (7) to have full period
2rw − 1 can be found in Matsumoto and Kurita [1992]. In the case where bothm and r are even, we have never encountered a full-period generator of type II,
despite making exhaustive searches for r = 4, 8, and 12. So we conjecture that one
cannot have full period in this case, but we have no proof.
4.4

Generators of type III

The generators of type III are based on the recurrence
vn = (I + Ha1 1 )vn−m1 ⊕ (I + Ha2 2 )vn−m2 ⊕ (I + Ha3 3 )vn−r(9)where Hi ∈ {L, R}, ai < w, and 1 ≤ i ≤ r. We denote this variant by (H1 , H2 , H3 ,m1 , m2 , a1 , a2 , a3 ). If (H1 , H2 , H3 , m1 , m2 , a1 , a2 , a3 ) provides a generator of maxACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·

9

imal period then so does (K1 , K2 , K3 , m1 , m2 , a1 , a2 , a3 ) where Ki 6= Hi and Ki ∈{L, R}, by the similarity matrix C = P.
5.

SPECIFIC GENERATORS PROPOSED BY MARSAGLIA

Marsaglia [2003] lists all parameters (a, b, c) that yield full-period type-I xorshift
generators with w = 32 and w = 64. We have checked his results and they arecorrect except for what seems to be a typo: for w = 32, (a, b, c) = (9, 5, 1) should
read (a, b, c) = (9, 5, 14). For the type-II and type-III generators, Marsaglia provides only a few sets of parameters. In the following subsections, we analyze theequidistribution and statistical properties of some generators he proposed.
5.1

Equidistribution properties

We computed the values of ∆1 for all type-I full-period generators with w = 32 andw = 64; they are listed in Panneton [2004]. These values range from ∆1 = 1 (goodequidistribution) to ∆1 = 153 (very bad). For example, a type-I generator qualified
by Marsaglia as one of his “favorites” is based on X2 with (a, b, c) = (5, 17, 13) andw = 32. This generator has ∆1 = 2, a good value. Another interesting example
is the type-I generator based on X4 with (a, b, c) = (7, 1, 9). This generator has∆1 = 1, the best value of all type-I generators with w = 32, but if we change X4
to X3 with the same parameters (a, b, c), i.e., simply change the order in which thetwo right xorshifts are performed, we get ∆1 = 56, which is the worst possible value
when w = 32. This illustrates the fact that the behavior of a generator based ona given triple (a, b, c) and given transition matrix Xi does not help predicting the
behavior of another generator with the same triple but a different choice of Xi .
All type-II and type-III generators proposed by Marsaglia with w = 32 haverather poor values of ∆1 . For example, the type-II generator ((I + R1 )(I + L2 ), (I +R4 ), 1) with r = 5 has ∆1 = 164, the type-II generator ((I + R1 )(I + L10 ), (I +R26 ), 1) with r = 3 has ∆1 = 81, and the type III generator (L, R, L, 1, 2, 6, 19, 3)
with r = 3 has ∆1 = 69.
5.2

Statistical Testing

We have applied the batteries of empirical statistical tests SmallCrush and Crush
implemented in the software package TestU01 [L’Ecuyer and Simard 2001b] toseveral xorshift generators proposed by Marsaglia [2003]. All tests included inthese batteries look for evidence against the null hypotheses H0 that the generator
produces i.i.d U (0, 1) random variables. These two batteries run in about one
minute and one hour, respectively, on a standard PC. None of the proposed xorshift
generators passed all these tests and most of them even failed the “baby” battery
SmallCrush in a spectacular way.
In what follows, we describe a few of these tests and give concrete illustrationsof the results.
The maximum-of-t test (see Knuth 1998, page 70) generates n sequences of tvalues in [0, 1] and computes the maximum X of the t values for each sequence.
The interval [0, 1] is partitioned into d segments in a way that under H0 , X falls inany given segment with probability n/d. The empirical and theoretical frequenciesof the d segments are then compared via a chi-square test statistic. The p-valueof this test is defined as the probability that the chi-square statistic takes a value
ACM Journal Name, Vol. V, No. N, Month 20YY.

10

·
Num.
1
2
3
4
5

F. Panneton and P. L’Ecuyer

Table II. List of tests from TestU01 used in this paper
Name
Parameters
Maximum-of-t
N = 10, n = 107 , τ = 0, d = 5 × 105 , t = 5
Birthday spacings
N = 5, n = 107 , τ = 14, d = 28 , t = 8, p = 1
Indep. Hamming weights
n = 108 , τ = 20, s = 10, L = 30
Matrix rankn = 106 , τ = 20, s = 10, L = 90
Matrix rankn = 50000, τ = 20, s = 10, L = 300

larger or equal to the one observed, under H0 . For a two-level test, this is repeated
N times and the empirical distribution of the N p-values thus obtained is comparedto the uniform distribution via a Kolmogorov-Smirnov test, whose p-value (at thesecond level) is taken as the final p-value of the test. In Crush, we find this testwith the parameters given in the first row of Table II (test number 1).
The birthday spacings test [Marsaglia 1985; L’Ecuyer and Simard 2001a] partitions [0, 1]t into dt subcubes of equal sizes, numbers them from 0 to k − 1 in anatural order, and throws n points into [0, 1]t using t successive uniform random
numbers for each point. Let I1 ≤ I2 ≤ · · · ≤ In be the (sorted) numbers of cells inwhich the points fall. The test computes the differences Ij+1 − Ij , for 1 ≤ j < n,and counts the number Y of collisions between these differences. Under H0 , Y isapproximately Poisson with mean λ = n3 /(4k). This process is repeated N times“independently” and the p-value of the test is defined as the probability that a
Poisson random variable with mean N λ takes a value larger or equal to the sum of
the N observed values of Y . In a slight generalization of the test, each output value
u of the generator (a real number between 0 and 1) is replaced by 2τ u mod 1 before
these numbers are used to determine the box number. The effect of this is to skip
the first τ bits of each output value of the generator and take the following ones.
This generalization applies to the other tests of Table II as well. Crush contains
this test with τ = 14, t = 8, d = 28 , n = 107 , and N = 5. So it uses bits 15 to 22
of 8 successive uniforms to determine the box number of each point, assuming thatbit 1 is the most significant bit. This is test number 2 in Table II.
The independence of Hamming weights test [L’Ecuyer and Simard 1999] takes s
successive bits, say bits τ +1 to τ +s, from each of 2n dL/se successive uniforms, andconcatenates these bits to construct 2n blocks of L bits. Let Xj be the Hamming
weight (the numbers of bits equal to 1) of the jth block, for j = 1, . . . , 2n. Eachvector (Xj , Xj+1 ) can take (L + 1) × (L + 1) possible values. The test counts
the number of occurrences of each possibility among the non-overlapping pairs{(X2j−1 , X2j ), 1 ≤ j ≤ n}, and compares these observations with the expected
numbers under H0 , via a chi-square test, after lumping together in a single classall classes for which the expected number is less than 10. This is test number 3 inTable II.
The matrix rank test generates n random L × L binary matrices, computes theirranks, and compares the empirical distribution of these ranks with their theoretical
distribution under H0 via a chi-square test [Marsaglia 1985; L’Ecuyer and Simard
2001a]. Each matrix is filled up line by line by taking s-bit blocks from L2 /ssuccessive uniforms (assuming that s divides L). Tests 4 and 5 in Table II are ofthis type. These tests skip the τ = 20 most significant bits of each output value and
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·

11

uses the s = 10 bits that follow. With a xorshift generator, the rank of the matrix
cannot exceed the degree k = rw of P (z), so we expect all xorshift generators to
fail this test when L > rw. But in many cases, we observed decisive failures for L
much smaller than kw.
We tested all full-period type-I generators with A1 = X1 , for w = 32 and w = 64.
All type-I generators with w = 32 and about half of those with w = 64 had a pvalue smaller than 10−300 for at least one test of SmallCrush. This can certainly
be called a spectacular failure! All but two of the generators of type II proposed
in Marsaglia [2003] failed a maximum-of-t test found in SmallCrush (with t = 6,
d = 105 , τ = 0, and n = 2 × 106 ) with a p-value smaller than 10−300 , and theother two failed other tests in SmallCrush or Crush. The generator qualified by
Marsaglia as one of his “favorites” had p-values smaller than 10−164 for all five testsof Table II.
The type II generator ((I + R1 )(I + L2 ), (I + R4 ), 1) with r = 5 has the p-valueof p = 3.6 × 10−9 in test 2. The one with ((I + R1 )(I + L10 ), (I + R26 ), 1) and r = 3
gets p = 1.9 × 10−165 for test 1 and p < 10−300 for test 3. The type-III generator
proposed by Marsaglia [2003], which has (L, R, L, 1, 2, 6, 19, 3) and r = 3, gets
p = 1.6 × 10−96 for test 1.
6.

A SEARCH FOR BETTER XORSHIFT GENERATORS

We made a search for good full-period generators of types II and III with respect to the equidistribution criterion ∆1 with the constraint that the number ofxorshifts cannot exceed 3, as for most generators proposed in Marsaglia’s paper.
We tested empirically the best generators we found, by applying SmallCrush andCrush [L’Ecuyer and Simard 2001b].
Table III lists the best generators of type II we found, with respect to ∆1 , inan exhaustive search among those with w = 32 and r = 2, 3, 4, 5, 8, 12, 25. The
best ones have reasonably good values of ∆1 if we compare with other well-known
generators. For example, for r = 25 and w = 32, we have one with ∆1 = 123,
which compares advantageously with the TT800 generator of Matsumoto and Kurita[1994], for which ∆1 = 261.
We applied SmallCrush and Crush to all generators listed in Table III. The
number p3 in the table is the p-value for test number 3 of Table II. The symbol means “smaller than 10−300 ” and a blank indicates a p-value between 0.01 and
0.99. For r ≤ 5, all but one generator fail this test spectacularly. This means that
there is significant dependence between the Hamming weights of successive numbers
produced by these generators. Many of these generators also failed the matrix rank
tests of Table II. For example, two of the generators with rw = 25 × 32 = 800
failed Test 5 and three of those with rw = 12 × 32 = 384 failed Test 4, both witha p-value smaller than 10−300 .
Three generators in Table III passed all the tests of SmallCrush and Crush: theseare generators 22, 28, and 29. In view of the weaknesses of the other generators of
the same family and same (or slightly smaller) period, we cannot recommend themfor simulation purposes.
We also made an exhaustive search for good generators of type III with respectto ∆1 , with w = 32 and r = 2, 3, 4, 5, 8, 12, 25, and tested the best ones with
ACM Journal Name, Vol. V, No. N, Month 20YY.

12·F. Panneton and P. L’Ecuyer

Table III.

Value of ∆1 for the best generators of type II with w = 32 and r = 2, 3, 4, 5, 8, 12, 25.

Num.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

r
2
2
2
2
3
3
3
3
4
4
4
4
4
5
5
5
5
8
8
8
8
12
12
12
12
25
25
25
25

m
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
3
3
3
3
1
1
5
5
1
5
9
9
7
2

H
(I + L11 )(I + L11 )(I + R8 )(I + L9 )(I + L11 )(I + R9 )(I + L23 )(I + L23 )(I + L23 )(I + L18 )(I + R7 )(I + L11 )(I + R7 )(I + L11 )(I + L17 )(I + L19 )(I + L11 )(I + R7 )(I + R7 )(I + L11 )(I + R6 )(I + L11 )(I + L9 )(I + R6 )(I + R6 )(I + L9 )(I + R13 )(I + L19 )(I + R14 )(I + L17 )(I + R7 )(I + L15 )(I + R11 )(I + L8 )(I + L21 )(I + R11 )(I + R6 )(I + L7 )(I + R8 )(I + L7 )(I + L7 )(I + R6 )(I + R8 )(I + L11 )(I + L11 )(I + R8 )(I + L20 )(I + L10 )G(I + R13 )(I + L19 )(I + L19 )(I + R13 )(I + L22 )(I + L17 )(I + R4 )(I + L13 )(I + L11 )(I + R7 )(I + R7 )(I + L11 )(I + R5 )(I + L13 )(I + L20 )(I + L19 )(I + R5 )(I + L12 )(I + R15 )(I + L7 )(I + L19 )(I + L20 )(I + L20 )(I + L20 )(I + L20 )(I + L8 )(I + L8 )(I + L10 )(I + L21 )(I + L6 )(I + L22 )(I + L18 )(I + L22 )(I + L18 )(I + L18 )(I + R13 )(I + L5 )(I + R13 )(I + L19 )

∆1
4
7
7
7
11
12
12
13
13
17
17
19
19
18
19
25
25
45
48
52
54
74
79
84
90
123
137
155
158

p3
10−6
10−9















SmallCrush and Crush. Almost all those with r ≤ 8 have failed at least one of
the tests, whereas all those with r ≥ 12 passed all the tests in these batteries.
The parameters of the latter generators are given in Table IV. All generators withr = 25 in that table have a smaller value of ∆1 than TT800, but slightly larger thanfor the best generators of type II. Although these generators behave better thanthose of type II in the statistical tests, we are not enthusiastic to recommend themfor simulation, because those of (already large) period lengths 28×32 − 1 = 2256 − 1
fail the tests decisively, which is not reassuring.
7.

INCREASING THE NUMBER OF XORSHIFTS

So far we restricted ourselves to generators with only three xorshift operations, as
in Marsaglia’s paper. An obvious idea for improving the statistical robustness is
to increase the number of xorshifts. In this context, it is important to note thatan a-bit left [right] xorshift does not modify the a least [most] significant bits. For
this reason, if we decide to have several xorshifts and if we care about the least
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators

Table IV.

·

13

Value of ∆1 for the best generators of type III with w = 32 and r = 12, 25.

Num.
31
32
33
34
35
36
37
38

r
12
12
12
12
25
25
25
25

m2
2
5
7
5
4
5
7
5

m1
3
11
9
10
10
24
24
16

Ha2 2
L7
L5
L18
L5
L21
L5
L5
L19

Ha1 1
R11
L18
R11
L18
R11
R11
R11
R11

Ha3 3
L21
R11
L5
R11
L7
L18
L18
L5

∆1
96
100
102
103
186
188
190
219

significant bits as well as the most significant ones, there should be a good balance
between the numbers of left and right xorshifts. On the other hand, our criterion∆1 gives more importance to the most significant bits than to the least significant
ones, because of the way equidistribution is defined. For this reason, computersearches based on ∆1 tend to return generators having mostly left xorshifts andvery few right xorshifts. To balance the two types of xorshifts in a search based on∆1 , we must impose constraints on their respective numbers. Another possibility
would be to modify the criterion ∆1 to take into account the equidistribution of
the low-order bits. We leave this as a topic for further investigation.
We performed a computer search for full-period xorshift generators with w = 32,
order r = 8, s = 7 xorshifts with at least 3 right ones, and the smallest possible ∆1 .
These generators have period length 2256 − 1 and the best we found have ∆1 = 9.
One of them has recurrence vn = (I + L9 )(I + L13 )vn−1 + (I + L7 )vn−4 + (I +R3 )vn−5 + (I + R10 )vn−7 + (I + L24 )(I + R7 )vn−8 and its characteristic polynomial
P (z) has 131 nonzero coefficients. A C implementation of this generator in given
in Figure 1.
In a similar search for s = 13 and the same values for the other parameters, thebest we found also has ∆1 = 9 and 129 nonzero coefficients in its characteristic
polynomial. Its recurrence is: vn = (I + L17 )vn−1 + (I + L10 )vn−2 + (I + R9 )(I +L17 )vn−4 + (I + R3 )vn−4 + (I + R12 )vn−5 + (I + R25 )vn−5 + (I + R3 )(I + R2 )vn−6 +(I + R27 )vn−7 + (I + R22 )vn−7 + (I + L24 )(I + R3 )vn−8 . These two generators pass
all the tests in Crush. On the other hand, they obviously fail matrix rank tests for
L > 256 and tests based on linear complexity, because their bit sequences follow
linear recurrences of order 256, modulo 2.
For comparison, we made a similar search without any constraint on the numberof right xorshifts, with s = 13 and r = 7. The best generator we found had ∆1 = 1
but a single right xorshift for 12 left xorshifts, and it failed test 4 in Table II.
Brent [2004b] recently proposed type-II xorshift generators based on the recurrencevi = (I + Rd )(I + Lc )vi−m ⊕ (I + Rb )(I + La )vi−r ,(10)with four xorshift operations, and the parameters given in Table V. We computed
the equidistribution of these generators and their values of ∆1 , given in the table, are
not particularly good. We submitted them to the Crush battery, and the first two
failed the matrix-rank tests of Table II, but the other ones (with r ≥ 8 and period
ACM Journal Name, Vol. V, No. N, Month 20YY.

14·

F. Panneton and P. L’Ecuyer

static unsigned int x[8];/* Generator’s state.*//* Initializes state.*/void initxorshift7 (unsigned int *init) {int j;for (j=0; j<8; j++) x[j] = init[j];}/* Advances by one step and returns a number in [0,1).*/
double xorshift7 (void) {static int k = 0;
unsigned int y, t;t = x[(k+7) & 0x7U];t = t ^ (t<<13);
y = t ^ (t<<9);t = x[(k+4) & 0x7U];y^= t ^ (t<<7);t = x[(k+3) & 0x7U];y^= t ^ (t>>3);
t = x[(k+1) & 0x7U];y^= t ^ (t>>10);t = x[k];t = t ^ (t>>7);y^= t ^ (t<<24);x[k] = y;k = (k+1) & 0x7U;
return ((double) y * 2.32830643653869628906e-10 );}Fig. 1.

A C implementation of a seven-xorshift generator

length ρ ≥ 2256 − 1) passed all the tests. Brent recognizes the potential weaknessesof these xorshift generators and proposes an implementation that incorporates a
combination with a Weyl generator, with the hope that this improves the quality.

Table V.

Values of ∆1 for the type II generators of the form (10) proposed by Brent.
r
m
a
b
c
d
∆1
2
1
17 14
12 19
7
4
3
15 14
12 17
34
8
3
18 13
14 15
58
16
1
17
15
13
14
142
32
15 19
11 13
16
141
64
59 19
12 14
15
465
128 95
17 12
13 15
845
132 67
15 14
13 18
1838
140 19
17 13
15 16
2038

To see how the number of xorshifts affects the speed of the generators, we implemented generators of types I, II, and III with three xorshifts, as well as thegenerators described above with 4, 7, and 13 xorshifts. To generate 109 (one billion) uniform random numbers in [0, 1) and add them up, it took approximately 32
seconds for the fastest generators (type I with three xorshifts) and 38 seconds withthe slowest one (13 xorshifts). So adding xorshifts does not slow down the generator
significantly in comparison with the time for the function call, transformation intoa real number, and adding the numbers. The timings with the Mersenne twisterof Matsumoto and Nishimura [1998] and the WELL generators of Panneton et al.
ACM Journal Name, Vol. V, No. N, Month 20YY.

Xorshift Generators·

15[2005] are about the same. These timings are for a 2.8Ghz Intel Pentium 4 processor
running Linux and the gcc compiler with the -O2 optimization flag.
8.

CONCLUSION

We have studied both theoretically and empirically the xorshift generators proposedby Marsaglia. These generators are fast, but not reliable, according to our analysis.
Those of type I are doomed right from the start because of their short period length.
The other ones also have problems when the number of xorshifts is restricted tothree. To get over these limitations, one may want to use a larger number ofxorshifts together with a long period. To get rid of the linear structure of the output
sequence and further improve the statistical robustness, these xorshift generators
could be combined with other RNGs from different classes. Preferably, this should
be supported by some theoretical analysis of the point set Ψt of the combined
generator, e.g., as in L’Ecuyer and Granger-Piché [2003].
ACKNOWLEDGMENTS

This work has been supported by the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC) Grant Number ODGP0110050, NATEQ-Québec grant
Number 02ER3218, and a Canada Research Chair to the second author. The firstauthor benefited from NSERC and NATEQ scholarships. We thank Richard Simard
for his help in improving the presentation and running the statistical tests.
REFERENCES
Brent, R. P. 2004a. Note on Marsaglia’s xorshift random number generators. Journal of Statistical Software 11, 5, 1–4. See http://www.jstatsoft.org/v11/i05/brent.pdf.
Brent, R. P. 2004b. Some uniform and normal random number generators. http://web.comlab.
ox.ac.uk/oucl/work/richard.brent/random.html.
Fushimi, M. and Tezuka, S. 1983. The k-distribution of generalized feedback shift register
pseudorandom numbers. Communications of the ACM 26, 7, 516–523.
Knuth, D. E. 1998. The Art of Computer Programming, Volume 2: Seminumerical Algorithms,
Third ed. Addison-Wesley, Reading, Mass.
L’Ecuyer, P. 1996. Maximally equidistributed combined Tausworthe generators. Mathematics
of Computation 65, 213, 203–213.
L’Ecuyer, P. 2004. Random number generation. In Handbook of Computational Statistics, J. E.
Gentle, W. Haerdle, and Y. Mori, Eds. Springer-Verlag, Berlin, 35–70. Chapter II.2.
L’Ecuyer, P. and Granger-Piché, J. 2003. Combined generators with components from different families. Mathematics and Computers in Simulation 62, 395–404.
L’Ecuyer, P. and Simard, R. 1999. Beware of linear congruential generators with multipliers of
the form a = ±2q ± 2r . ACM Transactions on Mathematical Software 25, 3, 367–374.
L’Ecuyer, P. and Simard, R. 2001a. On the performance of birthday spacings tests for certain
families of random number generators. Mathematics and Computers in Simulation 55, 1–3,
131–137.
L’Ecuyer, P. and Simard, R. 2001b. Testu01: Un logiciel pour appliquer des tests statistiquesà des générateurs de valeurs aléatoires. Unpublished software user’s guide.
Lidl, R. and Niederreiter, H. 1994. Introduction to Finite Fields and Their Applications,
Revised ed. Cambridge University Press, Cambridge.
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.
ACM Journal Name, Vol. V, No. N, Month 20YY.

16·

F. Panneton and P. L’Ecuyer

Marsaglia, G. 2003. Xorshift RNGs. Journal of Statistical Software 8, 14, 1–6. See http://www.jstatsoft.org/v08/i14/xorshift.pdf.
Matsumoto, M. and Kurita, Y. 1992. Twisted GFSR generators. ACM Transactions on
Modeling and Computer Simulation 2, 3, 179–194.
Matsumoto, M. and Kurita, Y. 1994. Twisted GFSR generators II. ACM Transactions on
Modeling and Computer Simulation 4, 3, 254–266.
Matsumoto, M. and Nishimura, T. 1998. Mersenne twister: A 623-dimensionally equidistributed
uniform pseudo-random number generator. ACM Transactions on Modeling and Computer
Simulation 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.
Niederreiter, H. 1995. The multiple-recursive matrix method for pseudorandom number generation. Finite Fields and their Applications 1, 3–30.
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., L’Ecuyer, P., and Matsumoto, M. 2005. Improved long-period generators basedon linear recurrences modulo 2. ACM Transactions on Mathematical Software. to appear.

Submitted November 23, 2004; Revised October 3, 2005.

ACM Journal Name, Vol. V, No. N, Month 20YY.