Elsevier

Computer Physics Communications

Volume 196, November 2015, Pages 161-165
Computer Physics Communications

The MIXMAX random number generator

https://doi.org/10.1016/j.cpc.2015.06.003Get rights and content

Abstract

In this paper, we study the randomness properties of unimodular matrix random number generators. Under well-known conditions, these discrete-time dynamical systems have the highly desirable K-mixing properties which guarantee high quality random numbers. It is found that some widely used random number generators have poor Kolmogorov entropy and consequently fail in empirical tests of randomness. These tests show that the lowest acceptable value of the Kolmogorov entropy is around 50. Next, we provide a solution to the problem of determining the maximal period of unimodular matrix generators of pseudo-random numbers. We formulate the necessary and sufficient condition to attain the maximum period and present a family of specific generators in the MIXMAX family with superior performance and excellent statistical properties. Finally, we construct three efficient algorithms for operations with the MIXMAX matrix which is a multi-dimensional generalization of the famous cat-map. First, allowing to compute the multiplication by the MIXMAX matrix with O(N) operations. Second, to recursively compute its characteristic polynomial with O(N2) operations, and third, to apply skips of large number of steps S to the sequence in O(N2 log(S)) operations.

Introduction

In  [1] it was proposed that k-mixing systems of Kolmogorov [2], [3], [4] may serve as a suitable random number generator. The particular system chosen was the one realizing linear automorphisms of the unit hypercube in RN:ui(t+1)=j=1NAijuj(t)mod1 where u[0,1). For the purposes of generating pseudo-random numbers with this method, one chooses the initial vector u(0), called the “seed”, with at least one non-zero component.

The entries of the matrix are integers: AijZ and subject to the following two conditions  [2], [3], [4], [1] on the defining matrix A:

  • detA=1

  • the eigenvalues λk of A must not lie on the unit circle, |λk|1 for all k=1,,N.

The first condition assures that the map defines a volume preserving automorphism. When the automorphism is viewed as defining a dynamical system with discrete time tZ, then the second condition assures that nearby trajectories diverge exponentially. There exists an everywhere dense, but discrete set of periodic trajectories all of which are unstable, and whose number as a function of the period τ asymptotically goes like ehτ/τ, where h is the Kolmogorov entropy of the system  [2], [5]: h=k:|λk|>1log|λk|. The auto-correlation decay time τ0 is related to the entropy as τ01/h. In this way, a connection is made between the chaotic dynamics of the system, its set of periodic trajectories and the entropy. In particular, the auto-correlation time is directly related to the empirical notion of randomness of the sequence. The second condition above states that none of the eigenvalues should lie precisely on the unit circle, and moreover the formula for the entropy demonstrates that it is desirable that most of the eigenvalues of the matrix should lie as far as possible away from it.

A particular matrix chosen in  [6] was defined for all N3: A=(11111112111113+s21111432111NN1N232). The MIXMAX matrix contains integer numbers and is defined recursively, since the matrix of size N+1 contains in it the matrix for size N. The only variable entry in the matrix is A32=3+s where s is some small “magic” integer, in many cases s=1 or s=0. For those N for which one or more eigenvalues lie on the unit circle for some s, one can choose another s.

The eigenvalues of the MIXMAX matrix are widely dispersed for all N, see Fig. 1. Thus, the spectrum of this system is multi-scale, with trajectories exhibiting exponential instabilities on all time-scales  [1].

Empirical evidence suggests that the largest eigenvalue appears to grow at least linearly with N, but we have been unable to obtain a strict bound on it. However, Kolmogorov’s entropy can be more conveniently calculated using the small eigenvalues as follows. Since the product of all of the eigenvalues is equal to the determinant, kλk=k|λk|=1,andklog|λk|=0, the entropy can be calculated equally well using the eigenvalues which are less than one by absolute value: h=k:|λk|<1log|λk| i.e. the entropy is equal to the area under the upper branch of the curve in Fig. 1 and also is equal to the area under the lower branch of the curve, taking into account that the vertical axis is already set to the logarithmic scale.

None of the eigenvalues of the matrix A are smaller by absolute value than 1/4, regardless of N, and a large number of them tend to cluster just above 1/4. Therefore the Kolmogorov entropy can be strictly bound from above by bounding the area under the lower branch of the curve as follows: h<Nlog|λmin|<Nlog(4).

We have also obtained a heuristic formula for the entropy which is more precise. The actual number of eigenvalues lesser than one by absolute value is asymptotically 2/3N, and logarithm of their value is observed to lie close to a parabola: log(λk)log(4)(1+(32N)2k2)for  k=1,,2N3. Adding these up, we get an approximate asymptotic formula, which also satisfies the strict bound above: Nlog(4)>h4/9Nlog(4). This estimate appears to be reasonable, i.e. for N=256 the actual value is h164.4 versus our estimate of h157.7.

Fig. 2 demonstrates graphically that some of the other popular generators do not satisfy the second requirement for randomness. In the case of RCARRY  [7] which is a slight modification of a Fibonacci-like recurrence modulo 224, this point has been made before by Lüscher  [8], and its failure was related to the weak mixing properties of its underlying matrix. Unfortunately, the Mersenne Twister (MT)  [9] has not been studied from this point of view. Our preliminary investigation indicates that the real eigenvalues of its characteristic polynomial are distributed very close to the unit circle, with the largest being less than |λmax|<1.0019, see Fig. 3. In this sense, the underlying dynamical system of the Mersenne Twister has one order of magnitude less entropy than RCARRY. This is related to the singular flaw acknowledged by the authors of the Mersenne Twister, which is that MT has a very long recovery time when it is seeded by a vector containing mostly zeros: the output is observed to be non-random even after outputting a million values. In this instance, the divergence of some trajectory is observed away from the origin. In fact, from the dynamical system point of view this is not merely a manifestation of an unlucky initialization, since any two nearby trajectories of the MT diverge very slowly and this is ultimately what causes the failure of MT in the statistical tests.

The total entropy of the MT system is approximately h4.8. On the basis of the known behaviour of RCARRY  [10], [11] and our investigation of the MIXMAX, it appears that an entropy of at least h50 is required for a generator to have a sufficiently random trajectory. In light of this, it is perhaps not surprising that the Mersenne Twister fails many tests in its pure form, and still fails some tests even when some additional tempering of the output is applied.

One can even conjecture that the sequence produced by the Mersenne Twister would equally benefit from the method of skipping proposed by Lüscher. However, the required amount of skipping would likely be prohibitively expensive. Moreover, it is clear that the recurrences based on primitive trinomials of even higher order such as those found in  [12] would exhibit even longer correlations.

Section snippets

Discrete case

In a typical computer implementation, the recursion (1) can be used to generate uniform random variables on the unit interval directly in floating point hardware. However, since actual floating point hardware has finite precision the sequence will tend to lie on a rational sublattice, with components of the vector xi being multiples of 2b where b is the number of bits available for the mantissa of the floating point unit, typically b=53 on current computers. There are at least two drawbacks to

Maximal period sequence

A key insight which allows to determine the maximal period of the sequence defined by (4) is that, in the case that the characteristic polynomial is primitive, Aq=p0Ifor  q=pN1p1, where as before p0 is the free term of the characteristic polynomial of A and is equal to the determinant of A. Therefore, the powers of the matrix A which are multiples of q are diagonal: Amq=(Aq)m=p0mI. Since the MIXMAX matrix has p0=1, therefore the maximum possible period is equal to q. We define two necessary

Efficient implementation

In this section, we give for completeness the formulae which allow the efficient implementation of the generator in actual computer hardware.

First, we present the formula which allows the efficient calculation of the recursion. Given the vector a with components ai, i=1,,N, a vector of partial sums b is formed according to b1=0,bi=bi1+aifori=2,,N. Then, the new vector is calculated: a1a1+bN,aiai1+bifor  i=2,,N. Finally, the correction due to the magic value is applied a3a3+sa2. It

Conclusion

The MIXMAX random number generator is currently made available by the author in a portable implementation in the C language at http://hepforge.org   [21]. We urge the reader to consult the README file included in the distribution at http://hepforge.org for details on the use of the generator. Also, an experimental implementation intended as part of a development version of the C++ ROOT library from CERN exists at  [22].

The generator outputs in double precision natively, with all 53 bits random.

Acknowledgements

I would like to thank Jan Ambjørn, Fred James and George Savvidy, the three people whose encouragement is chiefly responsible for my persistence over the years with implementing the generator in software and studying its theoretical properties. I thank J. Apostolakis, J. Harvey and L. Moneta at CERN for useful discussions.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 644121. We also

References (22)

  • G. Savvidy et al.

    On the Monte Carlo simulation of physical systems

    J. Comput. Phys.

    (1991)
  • M. Lüscher

    Comput. Phys. Comm.

    (1994)
    F. James

    Comput. Phys. Comm.

    (1994)
  • L.N. Shchur et al.

    The RANLUX generator: Resonances in a random walk test

    Internat. J. Modern Phys. C

    (1998)
  • F. James

    Finally, a theory of random number generation

  • A.N. Kolmogorov

    A new metrical invariant of transitive dynamical systems and automorphisms of Lebeg spaces

    Dokl. Acad. Nauk SSSR

    (1958)
  • V.A. Rokhlin

    On the endomorphisms of compact commutative groups

    Izv. Akad. Nauk

    (1949)

    On the entropy of automorphisms of compact commutative groups

    Teor. Ver. i Pril.

    (1961)
  • J.N. Franklin

    Deterministic simulation of random processes

    Math. Comp.

    (1963)

    Equidistribution of matrix-power residues modulo one

    Math. Comp.

    (1964)
  • D.V. Anosov

    Geodesic Flow on Closed Riemannian Manifolds of Negative Curvature

    (1967)
  • N. Akopov et al.

    Matrix generator of pseudorandom numbers

    J. Comput. Phys.

    (1991)
  • G. Marsaglia et al.

    Ann. Appl. Probab.

    (1991)
  • M. Matsumoto et al.

    Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator

    ACM Trans. Model. Comput. Simul.

    (1998)
  • Cited by (45)

    • Maximally chaotic dynamical systems

      2020, Annals of Physics
      Citation Excerpt :

      The examples of maximally chaotic systems were discovered and discussed in the earlier investigations by Artin, Hadamard, Hedlund, Hopf, Birkhoff and others [4–17], as well as in more recent investigations [18–31]. Here we shall introduce and discuss the classical- and quantum-mechanical properties of maximally chaotic dynamical systems, the application of the C-K theory to the investigation of the gauge [32–40] and gravitational systems [41–44], as well as their application in the Monte Carlo method [45–54]. The review is organised as follows.

    • A new Pseudo random number generator based on generalized Newton complex map with dynamic key

      2020, Journal of Information Security and Applications
      Citation Excerpt :

      One of the main elements of cryptography is the random numbers sequences which produced by certain generators. Random number generators (RNGs) are widely used in various fields of entertainment industry, video games, music production, economic models, artificial intelligence, testing digital communication systems, statistical simulation, software development, cryptography and so forth [6–8]. As an example of the applications of random number generators, we can point out Intel’s hardware implementation in 1999, which provided the first random number generator to add security to computer hardware [9].

    • HEJ 2: High energy resummation for hadron colliders

      2019, Computer Physics Communications
    View all citing articles on Scopus
    View full text