
Random Numbers and Deterministic Simulation
More recently, helping my colleague solve the question of the repeatability of a number of tests, I once again came across the task of simulating a device that generates sequences of random numbers. In this article I will talk about the difficulties that were discovered, which approach to resolving them was chosen, and how I checked that the new solution is better than the previous one. I note right away that the question of creating and verifying random sequences is very delicate, and almost no solution can be considered final in it. I would appreciate comments and corrections.
First, I will briefly talk about the subject area of hardware and pseudo-random generators, about their characteristics and requirements for them. Then I will turn to my particular case - the field of software simulation of computer systems, and to a specific problem that needed to be solved.

Sometimes the most reliable way to get a random number is to take it from a directory. Image Source: www.flickr.com/photos/nothingpersonal/337684768
Introduction: RNG, PRNG, RNG, DRNG
The task of obtaining a “random” sequence of numbers using a computer is surprisingly complex. Intuitively, this can be justified by the fact that the philosophy of digital systems and the architecture of modern computers are aimed at achieving absolute accuracy, at eliminating even a hint of uncertainty as a result of work. The stored and transmitted data are verified by checksums, electrical circuits are duplicated. Accident is expelled, so how can it be returned? And why?
In fact, many practical tasks require algorithms that use sufficiently long random sequences to solve. Simulation, genetic algorithms, numerical integration and solving systems of linear equations (Monte Carlo methods), computer games and modern cryptography. In each case, there is an understanding of which sequences can already be considered random and which are not yet. The used random number generators also differ (RNG, English RNG - random number generator).
It would seem that nature provides enough chaos, where can humanity keep up with it. It is enough to make a digitizing device for a certain process (thermal noise, flash in the sun, throwing a coin), connect it to a computer, and the thing is in the hat.
However, the practical use of “clean” hardware RNGs is limited by the following factors. Firstly, the speed of generating numbers may not be high enough, and applications can consume them so quickly that even modern interfaces such as PCI-Express will be loaded to the eyeballs. Secondly, the statistical characteristics of such a sequence may not be the best - the distribution is not right or changes over time, the width of the numbers is less than the required, the tendency to produce zeros more often than units (English bias), etc. Thirdly, the price of such a high-quality device can be quite high, which will limit its scope to areas with the most stringent requirements for randomness. Finally, the power consumption and size of the RNG may not be suitable for mobile devices.
For this reason, a combined approach is often used: the initial value (seed) is taken from a random source, and the entire sequence is generated from it using a nonrandom function:
(a 1 , s 1 ) = f (seed)
(a 2 , s 2 ) = f (s 1 )
(a 3 , s 3 ) = f (s 2 )
...
Here a 1 , a 2 , ... is the output pseudorandom sequence, and s 1 , s 2, ... is the internal state of the generator. The resulting construct is called a pseudo random number generator (PRNG - pseudo-random number generator).
The quality of the sequence obviously depends on the choice of the function f (), and sometimes also on the value of seed (for example, f (x) = a * x mod m gives a sequence of zeros at seed = 0). History knows many uses of bad features. I will give only two examples.
- a i + 1 = 65539 * a i mod 2 31 , seed - odd. This is the classic RANDU , widely used on UNIX from 1960-1970. Fails even the simplest tests (we will talk about them later).
- Donald Knuth [5] describes his attempt in 1959 to make a super-PRNG using an algorithm of 13 steps, one more confusing than the other. In practice, it turned out that the sequence immediately converged to one number, which was then transformed into itself.
Sources of randomness in computing platforms
Back to the hardware RNG. They have been present in the IBM PC architecture for some time and are offered both in Intel components and other vendors. A few examples below.
- At the end of the 20th century, a hardware RNG was built into the Intel 82802 Firmware Hub [3]. It worked properly, but there were complaints about it in terms of power consumption, generation speed, as well as compatibility with the ever-decreasing size of the process.
- For this reason, a new RNG was designed, located directly in the core of the CPU. It was presented in the Ivy Bridge microarchitecture processors as an RDRAND instruction [4]. Later it also began to be used for the RDSEED instruction. The differences between the two teams are written here . A couple of articles on Habré about RDRAND: habrahabr.ru/post/128666 habrahabr.ru/post/145188
- In modern platforms, its own RNG is located in the TPM module (English trusted platform module) and can be used both for the needs of cryptography and for other purposes. Details on the speed of operation and the quality of the issued numbers are described in [8].
- VIA also offers a hardware RNG as part of its VIA C3 processors [6].
General requirements for RNG
Since I work on software models of equipment, the problem arises of modeling the RNG present in them. Model i.e. ultimately, the function implementing the generator in C should satisfy the requirements put forward for physical RNGs. If this is not so, then the differences should be understood and explained by the features of the modeling process - such a question cannot be left to chance. In the end, every model is only an approximation of reality, it is important that it be good enough for practical needs, and its limitations be realized.
- An accident. The most vague criterion, and we will come back to it. But it’s already clear that it is wise to approach the choice of a function and test it with tests.
- Distribution. We assume that the required distribution is uniform U [0; 2 N -1], where N = 64. Modern generators produce exactly so many bits.
- Great period. Every sequence from the PRNG is looped due to the limited number of different internal states. This will not be noticeable if her period is very long. Our goal is a value of the order of 2 63 or higher.
- Generation rate. RDRAND can produce tens of millions of numbers per second. Simulation usually works much slower than real iron. We will not chase the high generation speed, however, you should not choose too complex functions if possible - you do not want to create a bottleneck.
- Predictability. This aspect is important for cryptographic applications. Our model is not used for such needs. In addition, cryptographically strong software PRNGs are usually very slow.
There is another requirement specific to simulation.
- Repeatability The state of the simulation at any time of its operation can be saved to disk at the save point (English checkpoint), transferred to another machine and restored there. In addition, it is necessary to support the possibility of reversing the passage of time , and the course of all processes, including those emerging from the PRNG numbers, must be repeated. This means that the simulation must store the value of the internal state of the generator.
Problems with rand ()
The very first and generally useless generator looked something like this:
// init seed int seed = ...; srand (seed); uint64 bad_rand () { return rand (); }
There is not even a hint of repeatability - although seed is clearly set at the beginning of the work, and the sequence will always be the same, when the simulation is rolled back to the previous state, it will not be repeated, since the state of the generator is uncontrollable during the work. In addition, bad_rand () has other problems similar to those described below. We will not return to it in the future.
Then the code for generating a 64 bit random number was rewritten like this:
uint64 unix_rand (uint64 * state) { srand (* state); uint32 r1 = rand (); srand (r1); uint32 r2 = rand (); * state = r2; uint64 result = ((uint64) r2 << 32) | r1; return result; }
Since the standard rand () from Libc returns an int, which in the worst case is 32 bits wide, two consecutive calls and a combination of them into one 64-bit number were used. Since rand () does not return the value of its internal state, we had to set it each time using srand ().
This approach was still bad for a number of reasons.
- Dependence of the upper half of the bits on the youngest. As the seed for the highest 32 bits, the lower ones act, i.e. parts of the number are connected by a functional dependence that can be detected, i.e. show randomness.
- RAND_MAX on my system was 2 31 -1. Those. actually rand () only produces 31 bits! And in the final number, only 62 bits change, the 31st and 63rd are always zeros.
- The most serious flaw is the uncertainty of the rand () function. In general, it cannot be said whether RAND_MAX will be the same on all systems. You can’t even claim that the implementation of rand () will be the same on Windows and Linux, will not depend on the bit depth of the OS, version of Libc, compiler or anything else. The same simulation running on different host systems will produce different results. A nightmare for debugging!
New solution
Finally, we brought the code to the following form:
const uint64 a1 = ...; const uint64 c1 = ...; const uint64 a2 = ...; const uint64 c2 = ...; uint64 two_lcg_rnd (uint64 * state1, uint64 * state2) { uint64 r1 = a1 * * state1 + c1; uint64 r2 = a2 * * state2 + c2; * state1 = r1; * state2 = r2; uint64 result = (r2 & 0xffffffff00000000ull) | (r1 >> 32); return result; }
- We used two linear congruent generators (English LCG) of the form f (x) = a * x + c mod 2 64 . The strengths and weaknesses of this type of LCG are well understood, and calculations in integers modulo a power of two are fast.
- The constants a1, a2, c1, c2 are chosen so that the LCG periods are large. Those interested can easily find their meanings on the Internet.
- In order for the generators to be independent, a1! = A2, c1! = C2, and two independent variables are used as the internal state. But these precautions do not guarantee independence; empirical verification is required.
- Each of the returned 64-bit LCGs uses only the most significant 32 bits. The lower bits of LCG have a significantly shorter period, so they are discarded.
Using our own function for the PRNG saves us from the most elusive error - non-determinism of simulation. But how random is its issue?
Empirical test
The principle of empirical analysis G (P) MF is the calculation for the input sequence of statistics that determines the truth of a hypothesis about the nature of its behavior. For example, you can check the equidistribution of pairs (triples, fours ...) of neighboring numbers, the absence of periods of monotonous growth and decline, etc. For a true RNG, none of these hypotheses should be confirmed with sufficient certainty. The more different statistics a particular PRNG passes, the higher its quality. If you are interested in comparing the quality of the RNG from various popular programs, I recommend reading [1].
There are several packages for RNG research, including Dieharder [2] and NIST SP 800-22 [7]. I used the TestU01 package [1]. It predefined three groups of tests that differ in the number of statistics checked in them: SmallCrush, Crush and BigCrush. Each of the subtests included in them gives the number p-value from the half-interval [0,1). To pass the subtest, it must be simultaneously both far from zero and from unity.
I used the medium-sized Crush set, which in all my experiments ran from half an hour to two hours. I launched the programs on such a system:
CPU: Intel Core (TM) i5-4300U CPU @ 1.90GHz
uname -a: CYGWIN_NT-6.3 hostname 1.7.31 (0.272 / 5/3) 2014-07-25 11:26 x86_64 Cygwin
For the second generator based on rand (), the results are deplorable: 139 out of 144 tests failed:
Hidden text
========= Summary results of Crush ========== Version: TestU01 1.2.3 Generator: unix_rand Number of statistics: 144 Total CPU time: 00: 46: 32.40 The following tests gave p-values outside [0.001, 0.9990]: (eps means a value <1.0e-300): (eps1 means a value <1.0e-15): Test p-value ---------------------------------------------- 1 SerialOver, t = 2 eps 2 SerialOver, t = 4 eps 3 CollisionOver, t = 2 eps 4 CollisionOver, t = 2 eps 5 CollisionOver, t = 4 eps 6 CollisionOver, t = 4 eps 7 CollisionOver, t = 8 eps 8 CollisionOver, t = 8 eps 9 CollisionOver, t = 20 eps 10 CollisionOver, t = 20 eps 11 BirthdaySpacings, t = 2 eps 12 BirthdaySpacings, t = 3 eps 13 BirthdaySpacings, t = 4 eps 14 BirthdaySpacings, t = 7 eps 15 BirthdaySpacings, t = 7 eps 16 BirthdaySpacings, t = 8 eps 17 BirthdaySpacings, t = 8 eps 18 ClosePairs NP, t = 2 3.2e-157 18 ClosePairs mNP, t = 2 3.2e-157 18 ClosePairs mNP1, t = 2 eps 18 ClosePairs mNP2, t = 2 eps 18 ClosePairs NJumps, t = 2 eps 19 ClosePairs NP, t = 3 3.2e-157 19 ClosePairs mNP, t = 3 3.2e-157 19 ClosePairs mNP1, t = 3 eps 19 ClosePairs mNP2, t = 3 eps 19 ClosePairs NJumps, t = 3 eps 19 ClosePairs mNP2S, t = 3 eps 20 ClosePairs NP, t = 7 1.8e-79 20 ClosePairs mNP, t = 7 1.8e-79 20 ClosePairs mNP1, t = 7 eps 20 ClosePairs mNP2, t = 7 eps 20 ClosePairs NJumps, t = 7 eps 20 ClosePairs mNP2S, t = 7 eps 21 ClosePairsBitMatch, t = 2 6.9e-65 22 ClosePairsBitMatch, t = 4 7.5e-143 23 SimpPoker, d = 16 eps 24 SimpPoker, d = 16 eps 25 SimpPoker, d = 64 eps 26 SimpPoker, d = 64 eps 27 CouponCollector, d = 4 eps 28 CouponCollector, d = 4 eps 29 CouponCollector, d = 16 eps 30 CouponCollector, d = 16 eps 31 Gap, r = 0 eps 32 Gap, r = 27 eps 33 Gap, r = 0 eps 34 Gap, r = 22 eps 35 Run of U01, r = 0 eps 36 Run of U01, r = 15 eps 37 Permutation, r = 0 eps 38 Permutation, r = 15 eps 39 CollisionPermut, r = 0 eps 40 CollisionPermut, r = 15 eps 41 MaxOft, t = 5 eps 41 MaxOft AD, t = 5 3.2e-157 42 MaxOft, t = 10 eps 42 MaxOft AD, t = 10 1.8e-79 43 MaxOft, t = 20 eps 43 MaxOft AD, t = 20 1 - eps1 44 MaxOft, t = 30 eps 44 MaxOft AD, t = 30 1 - eps1 45 SampleProd, t = 10 1 - eps1 46 SampleProd, t = 30 1 - eps1 47 SampleMean eps 48 SampleCorr 1 - eps1 49 AppearanceSpacings, r = 0 1 - eps1 50 AppearanceSpacings, r = 20 1 - eps1 51 WeightDistrib, r = 0 eps 52 WeightDistrib, r = 8 eps 53 WeightDistrib, r = 16 eps 54 WeightDistrib, r = 24 eps 55 SumCollector eps 56 MatrixRank, 60 x 60 eps 57 MatrixRank, 60 x 60 eps 58 MatrixRank, 300 x 300 eps 60 MatrixRank, 1200 x 1200 eps 62 Savir2 eps 63 GCD, r = 0 eps 64 GCD, r = 10 eps 65 RandomWalk1 H (L = 90) eps 65 RandomWalk1 M (L = 90) eps 65 RandomWalk1 J (L = 90) eps 65 RandomWalk1 R (L = 90) eps 65 RandomWalk1 C (L = 90) eps 66 RandomWalk1 H (L = 90) eps 66 RandomWalk1 M (L = 90) eps 66 RandomWalk1 J (L = 90) eps 66 RandomWalk1 R (L = 90) eps 66 RandomWalk1 C (L = 90) eps 67 RandomWalk1 H (L = 1000) eps 67 RandomWalk1 M (L = 1000) eps 67 RandomWalk1 J (L = 1000) eps 67 RandomWalk1 R (L = 1000) eps 67 RandomWalk1 C (L = 1000) eps 68 RandomWalk1 H (L = 1000) eps 68 RandomWalk1 M (L = 1000) eps 68 RandomWalk1 J (L = 1000) eps 68 RandomWalk1 R (L = 1000) eps 68 RandomWalk1 C (L = 1000) eps 69 RandomWalk1 H (L = 10000) eps 69 RandomWalk1 M (L = 10000) eps 69 RandomWalk1 J (L = 10000) eps 69 RandomWalk1 R (L = 10000) eps 69 RandomWalk1 C (L = 10000) eps 70 RandomWalk1 H (L = 10000) eps 70 RandomWalk1 M (L = 10000) eps 70 RandomWalk1 J (L = 10000) eps 70 RandomWalk1 R (L = 10000) eps 70 RandomWalk1 C (L = 10000) eps 71 LinearComp, r = 0 1 - eps1 71 LinearComp, r = 0 1 - eps1 73 LempelZiv 1 - eps1 74 Fourier3, r = 0 eps 75 Fourier3, r = 20 eps 76 LongestHeadRun, r = 0 eps 76 LongestHeadRun, r = 0 1 - eps1 77 LongestHeadRun, r = 20 eps 77 LongestHeadRun, r = 20 1 - eps1 78 PeriodsInStrings, r = 0 eps 79 PeriodsInStrings, r = 15 eps 80 HammingWeight2, r = 0 eps 82 HammingCorr, L = 30 eps 83 HammingCorr, L = 300 eps 84 HammingCorr, L = 1200 eps 85 HammingIndep, L = 30 eps 86 HammingIndep, L = 30 eps 87 HammingIndep, L = 300 eps 88 HammingIndep, L = 300 eps 89 HammingIndep, L = 1200 eps 90 HammingIndep, L = 1200 eps 91 Run of bits, r = 0 eps 91 Run of bits, r = 0 eps 92 Run of bits, r = 20 eps 92 Run of bits, r = 20 1 - eps1 93 AutoCor, d = 1 1 - eps1 94 AutoCor, d = 1 eps 95 AutoCor, d = 30 1 - eps1 96 AutoCor, d = 10 1 - 5.2e-10 ---------------------------------------------- All other tests were passed
For a new generator based on two independent LCGs, the picture is better:
Hidden text
========= Summary results of Crush ========== Version: TestU01 1.2.3 Generator: two_lcg_rnd Number of statistics: 144 Total CPU time: 00: 46: 13.21 The following tests gave p-values outside [0.001, 0.9990]: (eps means a value <1.0e-300): (eps1 means a value <1.0e-15): Test p-value ---------------------------------------------- 1 SerialOver, t = 2 eps 2 SerialOver, t = 4 eps 8 CollisionOver, t = 8 1 - 8.0e-13 13 BirthdaySpacings, t = 4 1.2e-82 15 BirthdaySpacings, t = 7 eps 16 BirthdaySpacings, t = 8 eps 17 BirthdaySpacings, t = 8 eps ---------------------------------------------- All other tests were passed
Only 7 tests out of 144 failed. Not bad, although for cryptographic needs I will not advise using this PRNG. I note that the individual "halves" of this generator also did not pass the same types of subtests.
Out of curiosity, I ran TestU01 (the longest BigCrush) on the sequence issued by the RDRAND instruction on the same CPU:
Hidden text
========= Summary results of BigCrush ========== Version: TestU01 1.2.3 Generator: rdrand Number of statistics: 160 Total CPU time: 15: 58: 06.92 The following tests gave p-values outside [0.001, 0.9990]: (eps means a value <1.0e-300): (eps1 means a value <1.0e-15): Test p-value ---------------------------------------------- 1 SerialOver, r = 0 eps 2 SerialOver, r = 22 eps 76 RandomWalk1 J (L = 1000, r = 0) 1.7e-4 ---------------------------------------------- All other tests were passed
Finding the answer to the question why even on a “truly” random generator TestU01 produces several failures, let us leave it for future researchers, although I have a hypothesis on this subject.
Conclusion
The theme of generating pseudorandom numbers and checking them for randomness is incredibly exciting. To all interested, I recommend reading, if you have not already done so, the corresponding chapter from the second volume of D. Knut [5]. The issue of using RNG for the needs of simulation, and not only computing systems, is discussed in [9].
Literature
- Pierre L'Ecuyer and Richard Simard. 2007. TestU01: AC library for empirical testing of random number generators. ACM Trans. Math. Softw. 33, 4, Article 22 (August 2007). DOI = 10.1145 / 1268776.1268777 simul.iro.umontreal.ca/testu01/tu01.html
- Robert G. Brown, Dirk Eddelbuettel, David Bauer. Dieharder: A Random Number Test Suite Version 3.31.1 www.phy.duke.edu/~rgb/General/dieharder.php
- Intel Corporation. Intel 810 Chipset Design Guide, June 1999. download.intel.com/design/chipsets/designex/29065701.pdf
- Intel Digital Random Number Generator (DRNG) Software Implementation Guide - software.intel.com/en-us/articles/intel-digital-random-number-generator-drng-software-implementation-guide
- Donald Knut. The art of programming. Volume 2. The resulting algorithms. Third Edition - 2011 - Williams. ISBN 978-5-8459-0081-4,
- Cryptography research, Inc. Evaluation summary: VIA C3 Nehemiah random number generator - 2003 - www.cryptography.com/public/pdf/VIA_rngsum.pdf
- National Institute of Standards and Technology. A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications - 2010 csrc.nist.gov/publications/nistpubs/800-22-rev1a/SP800-22rev1a.pdf
- Alin Suciu, Tudor Carean. Benchmarking the True Random Number Generator of TPM Chips. arxiv.org/ftp/arxiv/papers/1008/1008.2223.pdf
- Handbook of Simulation. Principles, Methodology, Advances, Applications, and Practice / Ed. J. Banks. - John Wiley & Sons, Inc., 1998. - ISBN 0-471-13403-1. - URL: books.google.com/books?id=dMZ1Zj3TBgAC