How good is this random number algorithm?

Question

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 2²⁶.

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.