# Algorithm for normal random numbers

## Abstract

We propose a simple algorithm for generating normally distributed pseudo random numbers. The algorithm simulates $`N`$ molecules that exchange energy among themselves following a simple stochastic rule. We prove that the system is ergodic, and that a Maxwell like distribution that may be used as a source of normally distributed random deviates follows in the $`N\mathrm{}`$ limit. The algorithm passes various performance tests, including Monte Carlo simulation of a finite 2D Ising model using Wolff’s algorithm. It only requires four simple lines of computer code, and is approximately ten times faster than the Box-Muller algorithm.

Pseudo random number (PRN) generation is a subject of considerable current interest. Deterministic algorithms lead to undesirable correlations, and some of them have been shown to give rise to erroneous results for random walk simulations , Monte Carlo (MC) calculations , and growth models . Most of the interest has been focused on PRN’s with uniform distributions. Less attention has been paid to non-uniform PRN generation.

Sequences of random numbers with Gaussian probability distribution functions (pdf’s) are needed to simulate on computers gaussian noise that is inherent to a wide variety of natural phenomena . Their usefulness transcends physics. For instance, numerical simulations of economic systems that make use of so called geometric Brownian models (in which noise is multiplicative) also need a source of normally distributed PRN’s . There are several algorithms available for PRN’s with Gaussian pdf’s. Some, such as Box-Muller’s algorithm, require an input of uniform PRN’s, and their output often suffers from the pitfalls of the latter . Robustness is therefore a relevant issue. In addition, Box-Muller’s algorithm is slow and can consequently consume significant fractions of computer simulation times . The comparison method demands several uniform PRN’s per normal PRN, and is therefore also slow . Use of tables  is not a very accurate method. Algorithms that are related, but not equivalent, to the one we propose here have been published , but they are somewhat cumbersome to use. In addition, no proof of their validity has been given.

We propose here a new algorithm for the generation of normally distributed PRN’s that is quite simple and fast. It is a stochastic caricature of a closed classical system of $`N`$ particles. Their velocities provide a source of PRN’s. We prove that, for any initial state, their pdf becomes Maxwellian in the $`N\mathrm{}`$ limit, after an infinite number of two-particle “collisions” take place. To this end, we first prove that our system is ergodic . The proof is not exceedingly difficult because our system is not deterministic. We also study its output as a function of $`N`$, and establish useful criteria for its implementation. Correlation test results are also reported.

For the motivation, consider numbers $`v_1,v_2\mathrm{}v_N`$, placed in $`N`$ computer registers, analogous to velocities of $`N`$ particles that make up a closed classical system in 1D. Pairs of registers $`ı`$ and $`ȷ`$, say, selected at random without bias, are to “interact” somehow, conserving quantity $`v_i^2+v_j^2`$. By analogy with the approach to equilibrium (i.e., to Maxwell’s velocities distribution) that is believed to take place in Statistical Physics, we expect that sufficient number of iterations will lead to an approximately Gaussian pdf of register values, from which the desired PRN’s may be drawn. (See also Ref. .) We define below the simplest interaction we can think of in order that (1) implementation on a computer be very fast, and (2) that we may be able to prove that a Gaussian pdf does indeed ensue.

Before the algorithm is implemented, all $`N`$ registers must be initialized to, say, $`v_ı=1`$ for all $`ı`$ satisfying $`1ıN`$, or all $`v_ı`$ may be read from a set of $`N`$ register values saved from a previous computer run, which we assume to fulfill $`v_ı^2=N`$. Let $`U(1,N)`$, $`U_ı(1,N)`$ be unbiased integer random variables, both in the interval $`[1,N]`$, except that $`U_ı`$ cannot equal $`ı`$. The algorithm follows:

$`ı=U(1,N);ȷ=U_i(1,N);`$ (1)   
$`v_ı(v_ı+v_ȷ)/\sqrt{2};`$ (2)   
$`v_ȷv_ı+\sqrt{2}v_ȷ`$ (3)   
The updated value of $`v_ı`$, from Eq. (2), is used in Eq. (3). After an initial warm up phase (see below), $`v_ı`$ and $`v_ȷ`$ may be drawn each time transformation (1-3) is applied. They are two independent PRN’s, each one with an approximately Gaussian pdf, with $`v_ı=0`$ and $`v_ı^2=1`$ for all $`ı`$, if $`N`$ is sufficiently large (see below). Transformation (1-3) may be thought of as a rotation of $`\pm \pi /4`$ with respect to a randomly chosen $`ıȷ`$ plane ($`+`$ and $``$ signs are for the two possible index orderings, $`ıȷ`$ and $`ȷı`$). Thus, quantity $`v_{\mathrm{}}^2`$ is conserved. Frequencies of events from sequences of $`10^6`$, $`10^8`$ and $`10^{10}`$ PRN’s generated with transformation (1-3), with $`N=1024`$, are exhibited in Fig.1.

We first explain why PRN’s genered by transformation (1-3) are expected to be normally distributed. Let $`𝐏_n(𝐯)`$ be the probability density at $`𝐯=(v_1,v_2,\mathrm{},v_N)`$, after transformation (1-3) has been applied $`n`$ times, on the $`(N1)`$-dimensional spherical surface $`𝒮_{N1}`$, of radius $`\sqrt{N}`$ given by $`N=_{\mathrm{}=1,\mathrm{},N}v_{\mathrm{}}^2`$. Let the single register pdf $`p(v)`$ be the $`n\mathrm{}`$ limit of $`p_n(v)`$, where $`p_n(v_1)=𝐏_n(𝐯)𝑑v_2𝑑v_3\mathrm{}𝑑v_N`$. We show further below that $`𝐏_n(𝐯)constant`$ over $`𝒮_{N1}`$ as $`n\mathrm{}`$. It then follows by integration that,

$$p(v)\left(1\frac{v^2}{N}\right)^{(N3)/2}.$$
(4)   
Clearly, $`p(v)C\mathrm{exp}(v^2/2)`$ in the $`N\mathrm{}`$ limit, which is the desired result.

We prove below, in three stages, that $`P_n(𝐯)`$ does indeed become homogeneous over spherical surface $`𝒮_{N1}`$, if $`N3`$, in the $`n\mathrm{}`$ limit. We first prove $`P_n(𝐯)P_n(𝐮)`$ as $`n\mathrm{}`$ if $`𝐯`$ and $`𝐮`$ are related. \[From here on, we say that points $`𝐯`$ and $`𝐮`$ are related if succesive transformations (1-3) of $`𝐯`$ can lead to $`𝐮`$.\] We then prove that the system’s “orbit” covers $`𝒮_{N1}`$ densely \[that is, that any point $`𝐯𝒮_{N1}`$ can be brought arbitrarily close to any other point $`𝐮𝒮_{N1}`$ by applying transformations (1-3) to $`𝐯`$ a sufficient number of times\]. Then, the desired result follows easily. It may help to place the significance of the proof that follows into proper perspective to note that if in Eq. (1) $`ȷU_i[1,N]`$ is replaced by $`ȷ=ı+1modN`$, the system becomes then non-ergodic, as can be easily checked numerically.

To start the proof, let kernel $`K(𝐯,𝐯^{})`$ be defined by $`P_{n+1}(𝐯)=K(𝐯,𝐯^{})P_n(𝐯^{})𝑑𝐯^{}`$, and let

$$F_n\{P_{n+1}^2(𝐯)P_n^2(𝐯)\}𝑑𝐯.$$
(5)   
Note first that $`F_n<0`$ implies that $`P_{n+1}(𝐯)`$ is more uniform than $`P_n(𝐯)`$, in the sense that $`𝑑𝐯[P_{n+1}(𝐯)\overline{P}]^2<𝑑𝐯[P_n(𝐯)\overline{P}]^2`$, where $`\overline{P}=1/𝑑v`$. It follows from the definition of $`K(𝐯,𝐯^{})`$ that

$$F_n=𝑑𝐯\{[𝑑𝐯_\mathrm{𝟏}K(𝐯,𝐯_1)P_n(𝐯_1)]^2P_n^2(𝐯)\}.$$
(6)   
Making use of the detailed balance condition, $`K(𝐯,𝐯^{})=K(𝐯^{},𝐯)`$, which our system satisfies, and the relation $`𝑑𝐯K(𝐯,𝐯^{})=1`$, Eq. (6) can be cast into,

$$F_n=\frac{1}{2}𝑑𝐯𝑑𝐯_\mathrm{𝟏}𝑑𝐯_\mathrm{𝟐}Q(𝐯,𝐯_1,𝐯_2),$$
(7)   
where, $`Q=K(𝐯,𝐯_1)K(𝐯,𝐯_2)[P_n(𝐯_1)P_n(𝐯_2)]^2`$. Therefore, in the $`n\mathrm{}`$ limit, $`P_n(𝐯)`$ becomes constant over each set in $`𝒮_{N1}`$ within which any two points $`𝐯,𝐮`$ are related.

We now prove that the system’s orbit covers $`𝒮_{N1}`$ densely. Let $`H_N`$ be the group of transformations in $`N`$ dimensions defined by Eqs. (1-3). We first show that any rotation in $`3`$D \[that is, any element of $`SO(3)`$\] can be approximated arbitrarily close by elements of $`H_3`$. The proof is extended to higher dimensions by induction. Note first that $`H_3`$ does not belong to the set of finite rotation groups in $`3`$D, and is therefore an infinite group. Let group $`SO(3)`$ be covered by spheres of radius $`ϵ/2`$ each. A finite number of them is sufficient, since the volume of $`SO(3)`$ is finite . It follows that there must be at least one sphere with two elements of $`H_3`$ in it, since $`H_3`$ has an infinite number of elements. Let these two elements be $`r`$ and $`s`$, and let $`g(𝐮,ϵ)`$ be element $`rs^1`$ of $`H_3`$, which is a rotation by angle $`ϵ`$ about some undetermined $`𝐮`$ axis. We will build elements of $`H_3`$ that are as near as desired to any given rotation. To this end, it is sufficient to show that it can be done for a set of infinitesimal generators of rotations . One such set is made up of infinitesimal rotations about three linearly independent axes. Consider axes $`𝐮_1`$, $`𝐮_2`$, and $`𝐮_3`$ that are obtained from $`𝐮`$ by rotations $`g(1,\pi /2)`$, $`g(2,\pi /2)`$, and $`g(3,\pi /2)`$ about each one of the coordinate axes by angle $`\pi /2`$. The correspondng infinitesimal rotations are given by , $`g(𝐮_ı,ϵ)=g(ı,\pi /2)g(𝐮,ϵ)g^1(ı,\pi /2).`$ This concludes the proof for 3 dimensions.

We now prove by induction that any element $`g(ıȷ,\alpha )`$, for the rotation about plane $`ıȷ`$, by angle $`\alpha `$ of the rotation group $`SO(N)`$ can be approximated as nearly as desired by an element $`g`$ of group $`H_N`$, for $`N>3`$. By hypothesis, any $`g(ıȷ,\alpha )`$, for $`ı,ȷ=1,2,\mathrm{},N`$ can be approximated by an element $`g`$ of $`H_N`$. We show now that $`g(ıN+1,\alpha )`$, for $`ı=1,2,\mathrm{},N`$, can also be approximated by elements of $`H_{N+1}`$. We take $`gH_N`$ within distance $`ϵ`$ of $`g_{ıȷ}(\alpha )`$. Now, since rotations preserve distances, it follows that $`g(ıN+1,\alpha )SO(N+1)`$, given by $`g(ıN+1,\alpha )=g(ıN+1,\pi /2)g(ıȷ,\alpha )g^1(ıN+1,\pi /2)`$ is within distance $`ϵ`$ of $`g^{}H_{N+1}`$, given by $`g^{}=g(ıN+1,\pi /2)gg^1(ıN+1,\pi /2)`$. This proves dense coverage in $`N3`$ dimensions. This is a stochastic generalization of Jacobi’s theorem  to more than two dimensions.

To conclude the proof that $`P_n(𝐯)constant`$ in the $`n\mathrm{}`$ limit, consider any two points $`𝐕`$ and $`𝐔^{}`$ as centers of disks $`𝒟_𝐕`$ and $`𝒟_𝐔^{}`$, both of radius $`r`$, in $`𝒮_{N1}`$. Since the system’s orbit covers $`𝒮_{N1}`$ densely for $`N3`$, it follows that a point $`𝐔`$ that is related to $`𝐕`$ exists arbitrarily close to $`𝐔^{}`$. Consider now disks $`𝒟_𝐕`$ and $`𝒟_𝐔`$. The fact that there exists at least one sequence of rotations in $`H_N`$ that take $`𝐕`$ into $`𝐔`$ implies that there exists at least one single rotation $`g`$ in $`H_N`$ that transforms $`𝐕`$ into $`𝐔`$. Since $`g`$ is a rotation, it transforms $`𝒟_𝐕`$ rigidly into $`𝒟_𝐔`$. It follows that $`𝑑𝐯P_N(𝐯)`$ over $`𝒟_𝐕`$ equals $`𝑑𝐮P_N(𝐮)`$ over $`𝒟_𝐔`$. Since $`r`$ is arbitrary, and $`𝐕`$ and $`𝐔^{}`$ are any two points in $`𝒮_{N1}`$, it follows that $`P(𝐯)`$ is constant over $`𝒮_{N1}`$ (except, perhaps, on a set of measure zero). This is the desired result. Ergodicity follows .

We next address the following practical issues: (1) how good an approximation to a Gaussian pdf of PRN’s is achieved with a necessarily finite set of $`N`$ registers; (2) how long must the warm up phase be.

It is convenient to rewrite Eq. (4) as follows,

$$p(v)e^{v^2/2}e^{g_N(v)/N}.$$
(8)   
where $`g_N(v)=v^2(3v^2/2)/2+𝒪(1/N)`$. $`N^1g_N(v)`$ is approximately the fractional deviation, $`\delta p(v)/p(v)`$, from Gaussian form if $`\delta p(v)/p(v)1`$. We have checked this behavior numerically. Clearly, the number of registers $`N`$ that must be used increases with the number $`M`$ of PRN’s one intends to generate. This is because the value of the largest PRN generated increases, on the average, with $`M`$. More precisely, the value of $`v`$ beyond which PRN’s are only generated with probability $`q`$ is approximately given by $`v^22\mathrm{ln}(M/vq)`$. Now, it follows from Eq. (8) that the fractional error $`\delta P/P`$ in the probability density at $`v`$ is approximately $`N^1v^2(3v^2/2)/2`$ for very large $`N`$. (It is pointless to require this error to be too small since a PRN is expected to be generated beyond $`x`$ with a small probability $`q`$.) It then follows that $`[\mathrm{ln}(M/qv)]^2N\delta P/P`$ must be satisfied by $`N`$. Thus, approximately $`10^4`$ registers are sufficient in order to generate as many as $`10^{15}`$ PRN’s, with a roughly $`10\%`$ error in the probability for the largest PRN in the sequence. For results obtained from a sequence of $`10^{10}`$ PRN’s generated with $`1024`$ registers, see Fig. 1.

Our algorithm must be applied a number $`n_pN`$ of times before it is ready for use unless all $`v_ı`$ are initialized with“equilibrium” values (stored from some previous computer run). The distribution of all register values then evolves towards equilibrium, as illustrated in Fig. 2. Deviations from equilibrium are statistically insignificant for $`n_p2`$ and $`N=1024`$, and for $`n_p4`$ and $`N=\mathrm{1\hspace{0.25em}048\hspace{0.25em}576}`$. Since $`n_p`$ is expected to increase as $`\mathrm{ln}N`$, $`n_p=8`$ should provide ample warm up for any forseeable applications.

The number of PRN’s that must be generated before each PRN in sequence $`v_1,v_2,\mathrm{},v_N`$ returns within distance $`r`$ from its initial value is exponential in $`N`$. More specifically, we estimate it to be $`(\tau /\sqrt{N})(1/r)^N`$ for $`N1`$, where $`\tau `$ is the period of the algorithm used to select $`ı`$ and $`ȷ`$ in Eq. (1). The estimation is based on $`P_n(𝐯)constant`$ over $`𝒮_{N1}`$ as $`n\mathrm{}`$. Thus, an effectively infinite recurrence time follows for any reasonable value of $`N`$.

Correlations between a finite number of PRN’s clearly vanish as $`N\mathrm{}`$, since $`ı`$ and $`ȷ`$ in Eq. (1) are supposedly independent PRN’s. We have searched for correlations in $`m`$ succesively generated PRN’s $`v_1,v_2,\mathrm{}v_m`$, for $`m=3,4,\mathrm{},6`$, performing a chi-square isotropy test over the corresponding $`m`$-dimensional space. An $`m`$-tuple $`𝐯=v_1,v_2,\mathrm{},v_m`$ was said to belong to the $`ith`$ cone, of $`1024`$ randomly oriented cones with axes $`𝐰_1,𝐰_2,\mathrm{},𝐰_{1024}`$, if $`0.99𝐯.𝐰_\mathbf{ı}1`$. No significant deviations from isotropy were observed for $`10^6`$ generated $`m`$-tuples.

Implementation of Wolff’s algorithm  in MC calculations of the Ising model’s critical behavior is a demanding test that some well known uniform PRN generators have failed . Large clusters are then flipped as a whole, and this tests correlations in very long sequences. We have used normal PRN’s generated by our algorithm as input into a MC simulation of an Ising system of $`16\times 16`$ spins at the critical temperature. \[For that, we note that $`v_ı^2+v_ȷ^2>2x`$ as often as $`u>\mathrm{exp}(x)`$ if $`v_ı`$ and $`v_ȷ`$ ($`u`$) are PRN’s with Gaussian (uniform) pdf’s, respectively.\] The energy obtained is shown in Fig. 3 as a function of the number of registers $`N`$. The following uniform PRN algorithms were used to select $`ı`$ and $`ȷ`$ in Eq. (1): ggl , R(250, 103,xor) , and RAN3 . We tried the latter two algorithms, which have been shown to lead by themselves to unacceptable results for the Ising model , in order to test our algorithm’s robustness. The results shown in Fig. 3 are gratifying.

Similarly, the specific heat $`c`$ and magnetization $`m`$ fluctuations data points obtained follow approximately the relations $`cc_0+8.4/N`$, and $`(\delta m)^2\chi _0+33/N`$, respectively, where $`c_0=1.497(1)`$ and $`\chi _0=0.5454(2)`$, in agreement with the known exact values.

Double precision is recommended. It prevents excesive drift of the sum $`v_ı^2`$ away from its assigned value. Even then, single precision accuracy is to be expected at the end of a sequence of some $`10^{16}`$ PRN’s, unless the sum is normalized several times during the run.

In summary, we have shown that implementation of Eqs. (1-3) provides a source of PRN’s with an approximately Gaussian pdf. Some $`10^4`$ registers (molecules) are sufficient for some purposes, but up to $`10^5`$ or more may be necessary for more demanding tasks. (Having to make a decision about the number of registers to be used may sometimes be an unwelcomed task. On the other hand, it is a virtue of the algorithm, that one can control, through the value of $`N`$, how close the output is to be from sequences of truly independent random numbers with Gaussian pdf’s.) Initial warm ups for arbitrary initial conditions are necessary; it is sufficient to let each register initially interact an average number of, say, 8 times. The system’s recurrence time was shown to be exponential in $`N`$, and therefore effectively infinite. Its behavior appears to be robust. The proposed algorithm runs an order of magnitude faster on computers than the most often used Box-Muller method . For a fortran code of our algorithm or other questions, please write JFF@Pipe.Unizar.Es.

Continuous help from Dr. Pedro Martínez with computer systems is deeply appreciated by JFF. We are indebted to Prof. P. Grassberger for an important suggestion. JFF and CC are grateful for partial financial support from DGICYT of Spain, through grants No. PB95-0797 and PB97-1080, respectively.