Fast K system generator of pseudorandom numbers

N.Z.Akopov,E.M.Madounts Yerevan Physics Institute,375036 Yerevan,Armenia

A.B.Nersesian Institute of Mathematics,375019 Yerevan,Armenia

G.K.Savvidy Physics Department,University of Crete,71409 Iraklion,Greece and

Institut für Theoretische Physik,D-6000 Frankfurt am Main, Germany

W.Greiner

Institut für Theoretische Physik, D-6000 Frankfurt am Main,Germany

Abstract

We suggest a new fast algorithm for the matrix generator of random numbers which has been earliar proposed in [1, 2]. This algorithm reduces N^{2} operation of the matrix generator to NlnN and essentially reduces the generation time. It also clarifies the algebraic structure of this type of K system generators.

In the articles [1, 2] the authors suggested the matrix generator of pseudorandom numbers based on Kolmogorov-Anosov K systems [3, 4]. This systems are the most stochastic systems, with nonzero entropy [3, 4, 5, 6, 7]. The properties of this new class of matrix generators were investigated with the criterion χ^{2} and the discrepancy D_{N} in different dimensions. In all cases it shows very good statistical properties [2] (for another examples see also [8]). Matrix generators based on different ideas are proposed in [9, 10].

In the recent article [11] this matrix generator was tested in the Monte- Carlo simulations of the Ising model in one and two-dimensions where we have defined control of the calculations. The same simulation was carried out with the RANECU generator. As it was shown in this particular case the quality of those generators are nearly the same, but the speed of the matrix generator, which we will call MIXMAX, is slower than that of the RANECU generator. It take place because the number of operations of the MIXMAX generator is of order N^{2}, where N is the dimension of the matrix.

The aim of this article is to suggest a new fast algorithm, which reduces N^{2} operations of MIXMAX to N ln N and allows to gain the generation time. It also clarifies the algebraic structure of this type of K system matrix generators.

The essence of the algorithm is that the almost T oeplitz structure of the matrix generator earliar proposed in [1, 2] rises up to circulant. The operations with the circulant matrix can be exchanged to the calculations with the diagonal matrix. The diagonalisation is performed by the direct and inverse discrete F ourier transformations and as it is well known, the fast F ourier transformation reduces N^{2} operations to N ln N. At this point the new algorithm accelerates the generation time.

Let us pass to the details of the algorithm. The matrix generator is defined as [1, 2],

X_{n+1} = A · X_{n}, (mod 1),  (1)

where A is N × N dimensional matrix [2]

A = ( \begin{matrix}2, 3, 4, ......., N , 1 \\ 1, 2, 3, ....., N − 1, 1 \\ 1, 1, 2, ....., N − 2, 1 \\ ................. \\ ................. \\ 1, 1, 1, ..., 2, 3, 4, 1 \\ 1, 1, 1, ..., 1, 2, 2, 1 \\ 1, 1, 1, ..., 1, 1, 2, 1 \\ 1, 1, 1, ..., 1, 1, 1, 1 \end{matrix} )  (2)

and X_{0} is an initial vector, whose components must be irrational. For example X_{0}^{(i)} = 1/\sqrt{π + i} and i = 1, .., N. The trajectory of the K system (1) X_{0}, X_{1}, X_{2}.... represents desired sequense of the random numbers [1]. This approach allows a large freedom in choosing of the matrices A for the K system generators and of the initial vectors [1].

Let’s present the matrix (2) in the form A = A_{1} − A_{2}, where

A_{1} = ( \begin{matrix}2, 3, 4, ..., N , N + 1 \\ 1, 2, 3, ..., N − 1, N \\ 1, 1, 2, .., N − 2, N − 1 \\ ................. \\ 1, 1, 1, ..., 1, 2, 3, 4 \\ 1, 1, 1, ..., 1, 1, 2, 3 \\ 1, 1, 1, ..., 1, 1, 1, 2 \end{matrix} ) , A_{2} = ( \begin{matrix}0, 0, 0, ..., 0 , N \\ 0, 0, 0, ..., 0, N − 1 \\ 0, 0, 0, ..., 0, N − 2 \\ ................. \\ 0, 0, 0, ..., 0, 0, 0, 4 \\ 0, 0, 0, ..., 0, 0, 1, 3 \\ 0, 0, 0, ..., 0, 0, 0, 2 \\ 0, 0, 0, ..., 0, 0, 0, 1 \end{matrix} ) .  (3)

So A_{1} is Toeplitz matrix and A_{2} is almost zero matrix. Let us extend the matrix A_{1} up to circulant Â_{1} of the dimension 2N × 2N

Â_{1} = ( \begin{matrix}2, 3, 4, ..., N , N + 1, 1 , 1 , .., 1 \\ 1, 2, 3, .., N − 1, N , N + 1, 1 , .., 1 \\ 1, 1, 2, .., N − 2, N − 1, N, N + 1, .., 1 \\ ................................. \\ ................................. \\ 5, 6, 7, ..., 1 , 1 , 1 , 1 , .., 4 \\ 4, 5, 6, ..., 1 , 1 , 1 , 1 , .., 3 \\ 3, 4, 5, .., N + 1, 1 , 1 , 1 , .., 2 \end{matrix} ) .  (4)

The following matrices

S_{1} = (E, 0) ,  S_{2} = ( \begin{matrix}E \\ 0 \end{matrix} ) ,  (5)

where E is identity matrix of the dimension N ×N, S_{1} is the matrix of the dimension N × 2N and S_{2} has dimension 2N × N allow to represent A_{1} in the form

A_{1} = S_{1}Â_{1}S_{2}.

Now the initial generator (1) can be written in the form,

X_{n+1} = (A_{1} X_{n} − A_{2} X_{n}), (mod 1),  (1a)

where

A_{1} X_{n} = S_{1}Â_{1}S_{2} X_{n}.  (6)

As it is known, the circulant matrix Â_{1} can be represented in the form:

Â_{1} = F^{−1}DF^{1},  (7)

where D is a diagonal matrix, F^{1} and F^{−1} are accordingly direct and inverse discrete Fourier transformations

F_{n,m} = (\frac{1}{\sqrt{2N}} exp \{ \frac{2πi}{2N}(n − 1)(m − 1) \})_{n,m=1}^{2N} ,  (8a)

and

D = diag (D_{1}.....D_{2N}) .  (8b)

Diagonal elements D_{1}, .., D_{2N} of the matrix (4) are counted by the formula:

D_{j} = 2 + 3r_{j} + ... + Nr_{j}^{N−2} + (N + 1)r_{j}^{N−1} + r_{j}^{N} + r_{j}^{N+1} + ... + r_{j}^{2N−1}, (9)

where r_{j} = exp( \frac{2πi}{2N}(j − 1) ), j = 1, .., 2N. Taking into consideration the last equation we may get the following formulae:

D_{1} = (N + 5)N/2,
D_{2} = −D_{1},
D_{j} = \frac{(2 + (−1)^{j−1}) (1 − (−1)^{j−1})}{(1 − r_{j})} + \frac{((−1)^{j−1}((N − 1)r_{j} − N) + r_{j})}{(r_{j} − 1)^{2}},  (10)

where j = 3, .., 2N. Now the algorithm of the matrix generator is realized by the scheme:

X_{n+1} = (S_{1} F^{−1} D F S_{2} X_{n} + A_{2} X_{n}) (mod 1)  (1b)

The representation (1b) of the MIXMAX (1) completely solves the problem.

This algebraic approach makes clear that at the ”heart” of the generator (1) [1, 2] is the T oeplitz matrix A_{1} in (3). Indeed one can check that the determinant of A_{1} is equal to one and the absolute values of all it‘s eigenvalues are different from one, therefore the matrix A_{1} also represents equally good the K system generator

X_{n+1} = A_{1} X_{n} = S_{1} F^{−1} D F S_{2} X_{n}

In fact the algebraic structure of this K system generators become more transparant in the sense that as we have seen one can embed the T oeplitz matrix A_{1} into the circulant of order 2N and decompose the last one into the sum of the powers of the basic permutation matrix,Ω

Â_{1} = 2 + 3Ω + 4Ω^{2} + ... + (N + 1)Ω^{N−1} + Ω^{N} + Ω^{N+1} + ... + Ω^{2N−1},  (9a)

where

Ω = ( \begin{matrix}0, 1, 0, ....., 0 \\ 0, 0, 1, ....., 0 \\ ............. \\ ............. \\ 0, 0, 0, ....., 1 \\ 1, 0, 0, ....., 0 \end{matrix} ) .  (11)

This also clearly demonstrates the fact that the growth of the matrix elements of the matrix A [1, 2] in (2) and A_{1} in (3) in the vertical direction produces reach spectrum of eigenvalues (9), the additional condition which is necessary to produce many − scale mixing of the system and to ensure a slower growth of the discrepancy D_{N}(A) [1].

It is also possible to modify the basic permutation matrix Ω in (11) to get a class of K system generators with very simple structure

A = ( \begin{matrix}0, 1 , 0 , ....., 0 \\ 0, 0 , 1 , ....., 0 \\ ............. \\ ............. \\ 0, 0 , 0 , ....., 1 \\ (−1)^{N+1}, a_{1}, a_{2}, .., a_{N−1} \end{matrix} ) .  (12)

In the last case it is easy to compute the characteristic polynomial of A

λ^{N} − a_{N−1} λ^{N−1} − ... − a_{1} λ + 1 = 0  (13)

and therefore it‘s eigenvalues λ_{1}, ..., λ_{N}

λ_{1} ⋯ λ_{N} = 1
............
λ_{1} + ... + λ_{N} = a_{N−1}.  (14)

This formulas allow to choose eigenvalues and then to construct matrix A for K system generators.

For example if N=4 and a_{1} = 0, a_{2} = 3, and a_{3} = 0, then

λ_{1} = \sqrt{\frac{3 + \sqrt{5}}{2}}, λ_{2} = −\sqrt{\frac{3 + \sqrt{5}}{2}},

λ_{3} = \sqrt{\frac{3 − \sqrt{5}}{2}}, λ_{4} = −\sqrt{\frac{3 − \sqrt{5}}{2}},  (15)

with an additional simplectic structure of A.

To convince that this algorithm produces the numbers of the same ”quality” as before [2] we check this algorithm for the matrix of the size 128 × 128. Bellow we present numerical results which were obtained for MIXMAX generator (1) realized though the algorithm (1b) together with RANECU generator. We have considered χ_{D}^{2} criterion when D = 1, .., 5

\begin{aligned}MIXMAX & D & χ^{2} \\ & 1 & 10.90 \\ & 2 & 88.89 \\ & 3 & 992.896 \\ & 4 & 9941.92 \\ & 5 & 100012.0 \end{aligned} \begin{aligned}RANECU & D & χ^{2} \\ & 1 & 7.625 \\ & 2 & 96.7 \\ & 3 & 990.52 \\ & 4 & 9976.64 \\ & 5 & 100416.25 \end{aligned}

The programm which compares generation time of the MIXMAX and RANECU shows that they work nearly with the same speed.

In conclusion the authors would like to thank F.James for interest and H.Niederreiter for provision of [9, 10]. This work was supported in part by the Alexander von Humboldt Foundation.

References

[1] G.K.Savvidy,N.G.Ter-Arutyunian-Savvidy, On the Monte-Carlo Simulation of Physical Systems, Preprint EPI-865(16)-86, Yerevan Jun.1986 [Russian]; J.Comput.Phys. 97, 566 (1991).

[2] N.Z.Akopov,G.K.Savvidy,N.G.Ter-Arutyunian-Savvidy, Matrix Generator of Pseudorandom Numbers, Preprint EPI-867(18)-86, Yerevan Jun.1986 [Russian]; J.Comput.Phys.97, 573 (1991).

[3] A.N.Kolmogorov, Dokl.Akad.Nauk SSSR 119, 861 (1958) [Russian].

[4] D.V.Anosov, Geodezicheskiye Potoki na Zamknutych Rimanovych Mnogoobraziyach Otrizatelnoi Krivizny (Geodesic flows on closed Riemannian manifolds of negative curvature) (Nauka, Moscow, 1967) [Russian].

[5] I.P.Kornfeeld, S.V.Fomin, Ya.G.Sinai, Ergodic theory (Springer-Verlag, New York, 1982) (Engl.transl).

[6] G.K.Savvidy, Nucl.Phys.B246, 302 (1984).

[7] G.K.Savvidy.Phys.Lett.130B, 303 (1983).

[8] R.O.Abramyan, N.Z.Akopov, G.K.Savvidy, N.G.Ter-Arutyunian-Savvidy, Sinai Billiards as a Pseudorandom Number Generator, Preprint EPI-922(73)- 86, Yerevan 1986 [Russian]; G.A.Galperin, N.I.Chernov.Billiardi i Chaos. Matematika i Kibernetika 5, 10 (1991), (Znanie, Moskva, 1991),[Russian].

[9] H.Niederreiter. Math. Japonica 31, 759 (1986).

[10] H.Grothe. Zufallszahlen und Simulation (Teubner, Stuttgard, 1986); Statist. Papers 28, 233 (1987).

[11] N.Z.Akopov, E.M.Madunts, G.K.Savvidy, A new matrix generator for lattice simulation, in Proceedings of Computing in High Energy Physics‘91 International Conference, pp.477-479 (Tsukuba, Japan, 1991).