# How good is this random number algorithm?

What follows are the random number generation routines from an ALGOL-60 computational math library, stored on a BESM-6 disk in a text form. The character encoding was with parity, imitating a punched tape character stream and with little to no line feeds; I've added line breaks and indentation:

``````_REAL _PROCEDURE RND(RANDOM);
_REAL RANDOM;
_BEGIN
_REAL A;
A:=SQRT(-2×LN(RANDOM));
RND:=A-(2.515517+A×(.802853+.010328×A))/
(1+A×(1.432788+A×(.189269+.001308×A)))
_END ;

_REAL _PROCEDURE RANDOM(Y);
_INTEGER Y;
_BEGIN
Y:=3125×Y;
Y:=Y-(Y÷67108864)×67108864;
RANDOM:=Y/67108864
_END ;
``````

The `RANDOM` function takes an integer and returns a real in [0; 1); the C equivalent would be ((3125*y) % 67108864) / 67108864.0, where 67108864 is 226.

The `RND` function takes a value (obviously positive, as its logarithm is taken) and produces the next supposedly random value. However, giving the resulting function to Wolframalpha shows that it is not particularly random, and its value range goes into the negative territory.

There were no comments or attributions. Does that function look like anything known in the field, but suffered transcription errors, or is it a complete garbage?

Update: The misunderstanding was due to the naming convention. The actual uniform number generator function is

``````_PROCEDURE UNR(X,A,V);
_INTEGER X;
_REAL A,V;
_BEGIN
_INTEGER I;
V:=0;
_FOR I:=1_STEP 1_UNTIL 5_DO _BEGIN
X:=3125×X;
X:=X-ENTIER(X/67108864)×67108864;
A:=X/33554432-1;
V:=V+A
_END ;
V:=V×0.774596;
V:=0.97×V×(1+0.01×V×V)
_END
``````

The same question about its provenance and quality stands.

• Wouldn't it make sense to ask it (as well?) at CS? It seems to be more of a CS/Math problem, doesn't it? Jan 14, 2022 at 4:15
• @Raffzahn If I don't get a full answer here in a few days, I'll ask there. Jan 14, 2022 at 6:33
• Are you just interested in the distribution, or in other aspects such as speed? "Good"ness depends greatly on what things you value... Jan 15, 2022 at 11:30
• UNR is not a uniform number generator but a standard normal number generator.
– WimC
Jan 15, 2022 at 14:04
• @TobySpeight How large a period is typically needed for the Monte-Carlo methods? A LFSR-like generator with a much larger period was previously known, and I would expect it to be faster because it does not use multiplication or division (the compiler was not optimizing), so I'm surprised. Jan 18, 2022 at 17:59

The short answer is that as a random number generator it's very poor algorithm, but on the other hand, the way it is used to generate random numbers with certain properties might be adequate for the task it was needed, so how good it is depends what you use it for.

The updated new function called UNR has something to do with generating uniformly distributed random numbers or something along those lines, but that's not important if you are simply asking about how the random number generarion algorithm works.

Generating the random numbers is a completely separate step from summing them up and generating the returned number with certain properties. The more random numbers are added together, the probabìlity density function of the result will get more like a Gaussian PDF. Apparently the amount of summed random numbers was reasonable balance between speed and quality of the generated number with the required distribution function.

The random numbers themselves are generated with the formula X(n+1) = (a•X(n) + c) mod m, so it is a type of linear congruential pseudorandom number generator (LCG).

Because c=0, it's a particular subtype called multiplicative congruential generator (MCG), sometimes called Lehmer RNG. The modulus m = 67108864 not even a prime number but a composite number that is a power of 2, which restricts the period to at most m/4, and the initial seed value for the generator must be carefully chosen or the period will be much less than that. It is also funny to see that the multiplier a = 3125 which is also not a prime.

Basically, MCGs are one of the simplest and earliest pseudorandom generators available, the Lehmer generator was published in 1951 and LCGs were published later in 1958. Their properties were known to be poor, but they were reasonably fast and simple as better alternatives were not yet available, as for example LFSRs were published in 1965.

• Note that the `UNR` procedure isn't a pure LCG/MCG: rather, it sums five successive outputs of the underlying MCG and then applies a cubic polynomial that (I assume) approximates the CDF of the resulting distribution in order to make the `V` output approximately uniform. (It does also return the last MCG output as an integer in `X`, and as a real scaled to the interval -1…+1 in `A`. I'm not sure if that's deliberate or just an artifact of the implementation.) Jan 14, 2022 at 12:38
• `UNR` is invoked, for example, in the `MONTECARLO` procedure. It is strange that the function used in the predecessor of the BESM-6 which looks faster and has a larger period, was not copied. Jan 14, 2022 at 16:30
• @LeoB. Why should they? Progress is a series of reinventing the wheel, so why searching for existing code to be reused when one can write in a few minutes? Jan 14, 2022 at 23:52

`RND` computes the inverse of the cumulative Gaussian distribution with mean 0 and standard deviation 1 (see, e.g., here). In other words, given a uniformly distributed random number in (0,1), it return a normally distributed random number.

The actual RNG must be elsewhere.

• Thank you; so the "bug" is in the naming convention. Jan 14, 2022 at 3:50
• I've found the actual RNG, please see update. Does it ring a bell? Jan 14, 2022 at 4:00

I guess that UNR stands for something like “unit normal” distribution, better known as the standard normal distribution. The algorithm gives a pretty decent sampling of that. Note that its range (for `V`) is restricted to the interval [-4.23, 4.23] but more than 99.998% of all expected values would be in that range for a standard normal variable.

As others have noted, the central loop adds five uniform variables in the range [-1, 1) to get a first approximation of a normal variable with a standard deviation of sqrt(5/3). The cubic is carefully chosen to further shape this closer to a standard normal distribution. It narrows the peak (bringing the variance close to one) and widens the tails. The first factor 0.774596 equals sqrt(3/5) to about 20 bits.

I could not readily find a reference, but I’d expect that this was a known method at that time. I’ll add more if I come across some other source.

Note 1. The number 3125 = 5^5 has order 2^24 (mod 2^26). This is the maximal order in the multiplicative group mod 2^26.

Note 2. This article states: “… in computations on the Russian computer BESM-6, for more than 20 years a sequence beginning with m0 = 1 while M = 2^40 and g = 5^17 has been used.” So a Lehmer generator (as used in the loop in this question) was common.

• Thank you for the link to the article; however, only one page is shown. What is the name of the publication and of the article? Jan 14, 2022 at 19:28
• @LeoB. Strange, I get all pages. It is ”A Pseudo-Random Number Generator for Personal Computers” by I.M. Sobol and Yu. L. Levitan. Computers and Mathematics with Applications 37 (1999) 33-40
– WimC
Jan 14, 2022 at 19:33
• Thank you; it is surprising that in a 1999 article the BESM-6 was still worth mentioning. Jan 14, 2022 at 22:53