IP/BBSR/93-52

Logistic Map: A Possible Random Number Generator

S. C. Phatak^{a} and S. Suresh Rao^{b}

Institute of Physics, Bhubaneswar-751005, India

Abstract

The logistic map is one of the simple systems exhibiting order to chaos transition. In this work we have investigated the possibility of using the logistic map in the chaotic regime (logmap) for a pseudo random number generator. To this end we have performed certain statistical tests on the series of numbers obtained from the logmap. We find that the logmap passes these tests satisfactorily and therefore it possesses many properties required of a pseudo random number generator.

PACS number(s): 05.45 +b, 05.20, 01.50.H, 02.70Lq

aphatak@iopb.ernet.in bsuresha@iopb.ernet.in

1 Introduction

A sequence of numbers that are chosen at random are useful in many different kinds of applications like simulation, sampling, numerical analysis, decision making, recreation etc. A sequence of truly random numbers is unpredictable and hence irreproducible. Such a sequence can only be generated by a physical process, for example, radioactive decay, thermal noise in electronic devices, cosmic ray arrival time etc. In practice, however, it is very difficult to construct physical generators which are fast enough and at the same time accurate and unbiased. Furthermore, one would like to repeat the calculation at will, for debugging or developing the program. Thus for most calculational purposes, pseudo random numbers (prn) have been introduced. Pseudo random numbers are numbers computed from a deterministic algorithm (hence it is called pseudo random or quasi-random) and therefore reproducible. Obviously these are not at all random in the mathematical sense, but are supposed to be indistinguishable from a sequence generated truly randomly. A good prn generator should possess long period, high speed and randomness.

Over the years various prn generators have been developed and can be broadly classified into the following categories [1]: those are, i.) Linear recurrence methods, ii.) Multiplicative congruential generators, iii.) Tausworthe generators and iv.) Combination generators. All these are bit based generators. These methods have various parameters or inputs and the period and the statistical properties of prn sequences sensitively depend on these parameters. The first two are known to have periods of about ∼ 10^{9} on a 32−bit machine. In the other two methods one can achieve much larger period (∼ 10^{170}) [1].

It is not easy to invent a foolproof source of random numbers. Generally a number of tests [2] are performed to test the randomness of the numbers generated by prn’s. In spite of these tests, one finds that a prn that has passed these tests may fail when applied to some physical applications: For example, in a recent study Ferrenberg et al. [3] have shown that even the high quality prns are biased under certain circumstances. Extensive Monte Carlo simulations by this group on an Ising model, for which exact answers are known, have shown that ostensibly high quality random number generators may lead to subtle, but dramatic, systematic errors for some algorithms, but not others. They traced the discrepancy to the correlations in the random numbers. Another recent study by Vattulainen et al. [4] find no such correlations. This only means that what is “random” enough for one application may not be random enough for another. The important criterion is that a specific algorithm must be tested together with the random number generator being used regardless of the tests which the generator has passed. This being the scenario, it may be useful to consider prn’s based on algorithms different from conventional algorithms.

In the present work we discuss a prn generator based on an inherently chaotic (random) algorithm and its statistical properties. Here we have employed the logistic map in the chaotic regime as a prn generator. We feel that such an effort is useful because it provides an entirely different method of producing prn’s. Also the algorithm used is very simple, so the generator is quite fast. There have been earlier attempts [5] to use the logistic map as a random number generator. For example, the logistic map in the chaotic regime was used to calculate the properties of a well-known random process; the invasion percolation problem. In this study, it was found that the static properties of percolating clusters, except for the percolation threshold, are correctly calculated but the dynamical properties are not. It is not clear from this work if this difference is due to non-uniform nature of the distribution or due to the correlation between successive numbers. The method we are adopting here differs from the earlier calculations in two respects. In the earlier works the random numbers were drawn from the logistic map which is known to have a distribution of ∼ \sqrt{x(1 − x)} and used in the simulation studies. We use a simple transformation to convert this distribution into a uniform distribution. Also, the earlier calculations make no attempt to either remove or study the effects of the correlations that exist between succesive numbers generated by the logistic map. We have done that. Here we draw the numbers from a uniform distribution (see below) and investigate the statistical properties. Our calculations show that the period of such a generator is certainly larger than 10^{9} (although theoritically infinite) if the computations are done in double precision. Our investigation of the distribution properties of this prn generator shows some peculiarities. These can be cured by introducing the τ–shift (to be explained below). In the following we shall first describe the logistic map and discuss the prn generator in Section 2. The results of different tests performed on prn’s thus generated are presented in Section 3 and the conclusions are given in Section 4.

2 The Logistic Map

Chaos in dynamical systems has been investigated over a long period of time [6]. With the advent of fast computers the numerical investigations on chaos have increased considerably over the last two decades and by now, a lot is known about chaotic systems. One of the simplest and most transparent system exhibiting order to chaos transition is the logistic map [7]. The logistic map is a discrete dynamical system defined by [7],

x_{i+1} = μx_{i}(1 − x_{i})  (2.1)

with 0 ≤ x_{i} ≤ 1. Thus, given an initial value (seed) x_{0}, the series x_{i} is computed. Here the subscript i plays the role of discrete time. The behavior of the series as a function of the parameter μ is interesting. A thorough investigation of logistic map has already been done[7]. Here, without going into detailed discussion, we simply note that,

• Eq.(2.1) has x = 0 and x = (μ − 1)/μ as fixed points. That is, if x_{i} = 0 or (μ − 1)/μ, then x_{i+1} = x_{i}.

• For μ < 1, x = 0 is an attractive (stable) fixed point. That is, for any value of the seed x_{0} between 0 and 1, x_{i} approaches 0 exponentially.

• For 1 ≤ μ ≤ 3, x = (μ − 1)/μ is an attractive fixed point.

• For 3 < μ < 4, the logistic map shows interesting behavior such as repeated period doubling, appearance of odd periods etc.

• For μ = 4 the logistic map is chaotic.

Since the chaotic behavior of the logistic map is of interest to us, we shall discuss the last point in detail. If we choose two seeds x_{0} and y_{0} = x_{0} + δx, with δx arbitrarily small, x_{i} and y_{i} differ by finite amount for large enough i. For example, if δx = 10^{−12}, x_{i} − y_{i} ∼ 0.1, for i ∼ 40. In fact, x_{i} − y_{i} grows exponentially with i and this is the definition of chaos. An analytic solution of logistic map exists for μ = 4. If we choose x_{i} = (1 − cos(θ_{i}))/2 and x_{i+1} = (1 − cos(2θ_{i}))/2 the logistic map for μ = 4 can then be defined by

θ_{i+1} = \{ \begin{matrix}2θ_{i}, & for θ_{i} < π/2 \\ 2π − 2θ_{i}, & for θ_{i} > π/2. \end{matrix}  (2.2)

Clearly, the map in terms of θ_{i} is given by stretching the line of length π to 2π and folding it. An examination Eq.(2.2) shows that one gets periodic series if the seed θ_{0} is a rational fraction of π. On the other hand, the series does not have periodicity if θ_{0} is an irrational fraction. Also, for any θ_{0} which is an irrational fraction of π, the set {θ_{i}} computed using Eq.(2.2) is distributed uniformly between 0 and π. It is this property of the logistic map that we intend to exploit for generating uniformly distributed random numbers and to study its statistical behavior.

Consider a set of numbers

y_{i} = \frac{cos^{−1} (1 − 2x_{i})}{π}  (2.3)

where,

x_{i+1} = 4x_{i}(1 − x_{i}).  (2.4)

From the discussion above, y_{i}’s are expected to be distributed uniformly between 0 and 1, except when x_{i} = \frac{1}{4}, \frac{1}{2}, \frac{3}{4} or 1. For these special values one gets periodic series. For any other rational fraction x_{0} we do not expect any periodicity. Thus it appears that, starting with a rational fraction x_{0} we can generate a set of uniformly distributed numbers {y_{i}} using Eqs(2.3,2.4). Furthermore since a small change in x_{0} produces large deviations in x_{i}’s (and therefore y_{i}’s) different initial values of x_{0} differing by small amount would produce different uncorrelated sets {y_{i}}. This seems to be a very good property for a random number generator to possess.

There are however two difficulties to be overcome before we can construct a random number generator from Eqs(2.3, 2.4). The first is due to truncation error in computers. We find that, when the calculations are done in single precision, x_{i} becomes 0 after about 5000 iterations. The exact value of i depends on x_{0} but it is ususally ∼ 5000. The reason for this is as follows. If x_{i} differs from 1/2 by a small amount ϵ (x_{i} = 0.5+ϵ), x_{i+1} = 1−ϵ^{2} and if ϵ^{2} < 10^{−7}, x_{i+1} is stored as 1 in the computer. So x_{i+n} = 0 for n ≥ 2. This difficulty is overcome by either doing the calculation in double precision (higher precision, if possible) or by modifying the algorithm suitably, when x_{i} is close to 0.5. We find that this problem does not arise for upto ∼ 10^{9} iterations when the calculation is done in double precision.

The second difficulty is that of correlation among the successive y_{i}’s. Eq(2.2) clearly shows that the succesive y_{i}’s are correlated, although a set of y_{i}’s are uniformly distributed in the interval [0,1]. The standard procedure of removing such a correlation is to shuffle the set of numbers obtained from logistic map. Our calculations show that the correlations have peculiar effect on the distribution properties and these persist even after shuffling (see later). On the other hand, the correlations between two numbers y_{i} and y_{i+τ} reduce as τ is increased from 1. We call these as τ-shifted numbers. The statistical tests have been performed on these τ-shifted set with τ ranging from 1 to 14. The following calculations have been done in double precision.

3 Distribution Tests

In order to test the prn generator (logmap) described in Section 2, we have performed certain tests, for various values of τ, described below. For comparison we have also done these tests on prn generator RAND– available on our machine and RANMAR– the generator proposed by Marsaglia and Zaman. These tests are:

1. Distribution test: here we verify Central Limit Theorem.

2. Moments Caclulation: We have calculated ⟨x^{n}⟩ and their variance σ_{n}.

3. χ^{2}-test: In the 1–d case, we have divided the interval [0, 1] into n equal bins and calculated χ^{2} and its distribution. For the 2–d case, we have divided the region [0, 1] ⊗ [0, 1] into n equal blocks and repeated the calculation of χ^{2} and its distribution.

In addition, we have also performed other tests [2] such as run-up test, n-tuple test etc. Results of these tests will not be presented here. The main reason is that the results of these tests are in concurrence with the above mentioned tests and it will not alter our conclusions. In the following, we report the results for logmap and RANMAR.

3.1 Distribution Test

Suppose that {x_{i}} is a sequence of mutually independent random variables that are governed by the probability density function P (x). Then the Central Limit Theorem asserts that, subject to certain conditions on the moments of P (x) [8] the variable y_{N} = ∑_{i=1}^{N} x_{i}, in the limit of large N, are distributed normally: i.e.,

P_{N}(y) = \frac{1}{\sqrt{2πσ^{2}}} exp −\frac{(y − μ)^{2}}{2σ^{2}} ,

where μ = N⟨x⟩ and σ^{2} = N(⟨x^{2}⟩ − ⟨x⟩^{2}).

We have obtained the distribution function of the numbers obtained from logmap for N = 24 for various τ–shifts. In fig.(1), the soild line indicates the exact distribution and diamond dots indicate the distribution obtained from the logmap with τ = 1, 2 and 3 calculation. From fig.(1), one can see that for τ = 1 (sequential case), the distribution from the logmap is skewed a little. Though the distribution for τ = 2 is better agreed than the sequential case, still it is far from satisfactory. For τ = 3 and greater, the distribution is in excellent agreement with the theoretical one. Thus this test shows that, although the successive numbers are correlated, a set of numbers obtained by picking every alternate or every third number, would behave as a set of random numbers [9].

3.2 Moments

For a given set of τ–shifted N numbers generated from logmap, we calculate n^{th} moment ⟨x_{τ}^{n}⟩ as

⟨x_{τ}^{n}⟩ = \frac{1}{N} ∑_{i=1}^{N} x_{i+τ}^{n}.  (3.1)

We then repeat this test for M such sets and caclulate the global average \overline{⟨x_{τ}^{n}⟩} and variance σ_{n} as

\overline{⟨x_{τ}^{n}⟩} = \frac{1}{M} ∑_{j=1}^{M} ⟨x_{τ}^{n}⟩_{j}  (3.2)
σ_{n}(τ) = \frac{1}{M} ∑_{j=1}^{M} ⟨x_{τ}^{n}⟩_{j}^{2} − \overline{⟨x_{τ}^{n}⟩}^{2}  (3.3)

We have varied τ from 1 to 14 for two different values of N–1024 and 8096 and M is held fixed at 4000. For random numbers uniformly distributed between 0 and 1 \overline{⟨x^{n}⟩}^{th} = \frac{1}{n+1} and σ_{n}^{th} = \frac{n^{2}}{(2n+1)(n+1)^{2}} . The results for the prn generator based on logistic map are shown in Table 1.

An inspection of Table 1 shows that \overline{⟨x_{τ}^{n}⟩} agrees well with the theoretical values for different n’s and τ’s. The agreement is comparable with the results of RANMAR. But σ_{n}(τ) departs significantly from the theoretical value for τ = 1, and the departure systematically increases with the increase in n. The variance is calculated by subtracting two almost equal numbers, so some loss of accuracy is expected. But the departure is much larger. On the other hand, for τ larger than 4, the calculated variances are close to the theoretical values and are comparable with those obtained from RANMAR.

3.3 χ^{2}-Test

In 1–D χ^{2} test, the interval [0,1] is divided into n equal parts and the numbers r_{i} falling in i^{th} interval, out of a set of N numbers is calculated. For a uniform distribution, χ^{2} is defined by,

χ_{n}^{2} = ∑_{i=1}^{n} \frac{(r_{i} − N/n)^{2}}{N/n}  (3.4)

The test is repeated for M such sets and the distribution of χ_{n}^{2} is obtained. Theoretically, the χ^{2} thus obtained should have a χ^{2} distribution for n − 1 degrees of freedom:

P_{th}(χ_{n}^{2}) = \frac{(χ_{n}^{2})^{\frac{n−3}{2}} e^{−χ_{n}^{2}}}{∫ dχ_{n}^{2}(χ_{n}^{2})^{\frac{n−3}{2}} e^{−χ_{n}^{2}}} .  (3.5)

The calculation is done for N = 10240, M = 4000, n is varied from 2 to 256 and τ is varied from 1 to 14. The results for n = 4, 64 and 256 are displayed in figs(2a,2b and 2c). The corresponding χ_{n}^{2} for RANMAR is also given for comparison.

Consider τ = 1 or the sequential case first. The calculated χ^{2} distribution differs significantly from the theoretical one. This is somewhat surprising, since the set of numbers obtained from the logistic map are are uniformly distributed. We have confirmed that shuffling does not mitigate this problem. Thus the correlations between successive numbers are probably responsible for this behavior.

As τ is increased, the agreement between calculated and theoretical distributions improves. From fig.(2a), we find that for n = 4, τ = 2 distribution already agrees reasonably with the theoretical distribution. On the other hand, for n = 256, one gets a good agreement for τ ∼ 6 or larger. This clearly shows that, in order to have ‘randomness’ on a finer scale (corresponding to smaller bin size, or larger n), the τ-shift must be larger.

For 2–D χ^{2}–test, we choose a pair of succesive numbers (x and y coordinates) and determine how they are distributed in n equal-area blocks covering a square of unit side. The χ^{2} distributions are calculated as discussed in the 1–D case and the results for n = 4, 64 and 256 are presented in figs(3a, 3b and 3c). These results follow the same pattern as that of 1–D case with one difference. For τ = 1, the calculated χ^{2} distribution is nowhere near the theoretical distribution for n > 4. This is simply because the succesive numbers are highly correlated. This can be explained as follows. For the sequential case, when the unit square is divided into four equal blocks, each block would have some number of points. Hence for n = 4, one does not expect any abnormal behavior in the χ^{2} distribution. However, this is not the situation when the unit area is divided into more than 2 × 2 blocks; it would so happen that some of the blocks do not contain any points at all. The contribution of these blocks to the χ^{2} (see eq.(3.4)) would be N/n, thus pushing the value of χ^{2} higher. However, with a τ–shift of 4 or larger, these correlations are more or less wiped out and one gets reasonably good agreement between calculated and theoretical distributions.

4 Summary and Conclusions

Various tests on the series of numbers obtained from the logmap have been performed in this work. These tests bring out certain peculiarities which have not been noted before. We notice that,

1. the distribution test does not satisfy the central limit theorem for τ = 1. The correlations between succesive numbers, which are known to exist, are probably responsible for this. But for τ > 2, the agreement with the central limit theorem is excellent.

2. The moments calculation comfirms the result of the distribution test. That is the calculated moments agree with their theoretical values, even for τ = 1. The variances of these moments disagree with their theoretical values for τ = 1. However, for τ > 4, the moments as well as variances agree with the theoretical values. This shows that the correlations between succesive numbers play a subtle role in moments calculation and removal of these correlations (by τ–shifting) is essential.

3. the χ^{2}–tests show that, for τ = 1 the χ^{2}–distributions systematically differ from the theoretical distribution. The agreement between calculated and theoretical distributions is improved by increasing τ. Thus, one must use τ–shifted numbers (with larger value of τ for smaller interval-size) to obtain acceptable χ^{2}–distribution.

As far, computer time, the logmap is about three times slower than the RANMAR However, in the calculations where prn’s are used, a relatively small fraction of time is spent in generating random numbers. Therefore the relative slowness of the logmap is not a big handicap. Furthermore, independent sets of random prn’s can be easily obtained from the logmap by choosing different values of seeds x_{0}’s. This is so because two seeds differing by a small amount would generate entirely different series after a few iterations. This property makes the logmap easily adaptable to parallel processing.

To conclude, we note that τ–shifted logmap satisfies some of the elementary tests a pseudorandom number generator must pass. The logmap being based on a physically chaotic process, the calculations where randomness, as opposed to computer time is important, it is advantageous to use logmap.

Acknowledgements

We gratefully acknowledge many useful discussions with Profs. R. E. Amritkar, S. M. Bhattacharjee, A. M. Jayannavar and A. Khare.

References

[1] For a recent review on random number generators see e.g., James, F. Comput. Phys. Commun. 60, 329 (1990); G. Marsaglia, B. Narasimhan and A. Zaman, ibid., 345 (1990).

[2] D. E. Knuth, The Art of Computer Programming, vol. 2: Seminumerical Algorithms 2nd ed. (Addison Wesley, Reading, 1981).

[3] A. M. Ferrenberg, D. P. Landau and Y. J. Wong, Phys. rev. Lett. 69, 3382 (1992).

[4] I. Vattulainen, K. Kankaala, J. Saarinen and T. Ala-Nissila, 1993, Helsinki University preprint, HU-TFT-01-93.

[5] S. Prakash, C. K. Peng and P. Alstrøm, Phys. Rev. A43, 6564 (1991); C. K. Peng, S. Prakash, H. J. Herrmann and H. E. Stanley ibid. A42, 4537 (1990).

[6] E. N. Lorentz, J. Atmos. Sci. 20, 130 (1993); A. Khare, Chaos, to appear in Physics Education (India); M. Taboor, Chaos and Integrability in Nonlinear Dynamics (John Wiley and Sons, N. Y.)

[7] Feigenbaum M. J., J. Stat. Phys. 19, 25 (1978).

[8] See e.g., the appendix of A. I. Khinchin, Mathematical Foundations of Statistical Mechanics, (Dover,1949)– translated from the Russian edition by G. Gamow, for details.

[9] We find that for τ = 1 and with shuffling, the agreement with the theoretical distribution is very good. This indicates that the departure for τ = 1 case is due to the correlation between successive numbers. However, the shuffling does not improve the results for moments and χ^{2} tests.

Figure Captions

fig.(1) Central Limit Theorem Verification for 24 numbers. The solid curve is obtained from the exact distribution P_{N}(y) = \frac{1}{\sqrt{2πNσ^{2}}} exp −\frac{(y−Nμ)^{2}}{2Nσ^{2}} , where N = 24, μ = 1/2 and σ^{2} = 1/12. The diamonds with dot are obtained from the logmap with τ = 1, 2 and 3 (as indicated by the label).

fig.(2) The χ^{2} distribution of the logmap for the 1−dimension case. The smooth curve is obtained from the exact distribution, eqn.(3.5) and the histogram is obtained from the logmap and the RANMAR, labelled appropriately. a) for 4 bins, b) for 64 bins and c) for 256 bins.

fig.(3) The χ^{2} distribution of the logmap for the 2−dimension case. The smooth curve is obtained from the exact distribution, eqn.(3.5) and the histogram is obtained from the logmap and the RANMAR, labelled appropriately. a) for 4 bins, b) for 64 bins and c) for 256 bins.

n  τ = 1 τ = 3 τ = 6 RANMAR exact
1 ⟨x⟩ .50010 .50006 .50000 .49993 .50000
σ_{1} .08400 .08460 .08626 .08155 .08333
2 ⟨x^{1}⟩ .33340 .33341 .33332 .33324 .33333
σ_{2} .03736 .08784 .09102 .08736 .08889
3 ⟨x^{3}⟩ .25004 .25007 .24999 .24991 .25000
σ_{3} .01392 .07828 .08111 .07964 .08036
4 ⟨x^{4}⟩ .20002 .20006 .20000 .19992 .20000
σ_{4} .00518 .06839 .07092 .07086 .07111
5 ⟨x^{5}⟩ .16667 .16672 .16667 .16660 .16667
σ_{5} .00273 .06002 .06240 .06305 .06313
6 ⟨x^{6}⟩ .14285 .14291 .14286 .14280 .14286
σ_{6} .00276 .05318 .05550 .05642 .05651
7 ⟨x^{7}⟩ .12498 .12505 .12500 .12495 .12500
σ_{7} .00363 .04759 .04990 .05087 .05104
8 ⟨x^{8}⟩ .11109 .11115 .11112 .11107 .11111
σ_{8} .00469 .04298 .04528 .04620 .04648
9 ⟨x^{9}⟩ .09997 .10004 .10001 .09997 .10000
σ_{9} .00567 .03913 .04143 .04225 .04263

Table 1: ⟨x^{n}⟩ and σ_{n} for logmap with τ = 1, 3 and 6 and RANMAR for n = 1, 2, ⋯, 10. N = 8196 and M = 4000.