11

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.

5
  • 1
    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?
    – Raffzahn
    Jan 14 at 4:15
  • 1
    @Raffzahn If I don't get a full answer here in a few days, I'll ask there.
    – Leo B.
    Jan 14 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 at 11:30
  • UNR is not a uniform number generator but a standard normal number generator.
    – WimC
    Jan 15 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.
    – Leo B.
    Jan 18 at 17:59

3 Answers 3

14

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.

3
  • 3
    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 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.
    – Leo B.
    Jan 14 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?
    – Raffzahn
    Jan 14 at 23:52
3

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.

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

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.

3
  • 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?
    – Leo B.
    Jan 14 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 at 19:33
  • Thank you; it is surprising that in a 1999 article the BESM-6 was still worth mentioning.
    – Leo B.
    Jan 14 at 22:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.