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.
    1. 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).
    2. 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.
    1. 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.
    2. 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.
    3. 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


    1. 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
    2. 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
    3. Intel Corporation. Intel 810 Chipset Design Guide, June 1999. download.intel.com/design/chipsets/designex/29065701.pdf
    4. Intel Digital Random Number Generator (DRNG) Software Implementation Guide - software.intel.com/en-us/articles/intel-digital-random-number-generator-drng-software-implementation-guide
    5. Donald Knut. The art of programming. Volume 2. The resulting algorithms. Third Edition - 2011 - Williams. ISBN 978-5-8459-0081-4,
    6. Cryptography research, Inc. Evaluation summary: VIA C3 Nehemiah random number generator - 2003 - www.cryptography.com/public/pdf/VIA_rngsum.pdf
    7. 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
    8. Alin Suciu, Tudor Carean. Benchmarking the True Random Number Generator of TPM Chips. arxiv.org/ftp/arxiv/papers/1008/1008.2223.pdf
    9. 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





    Also popular now: