Pseudo-random number generation algorithm in Redis source code

Guide----------------------------------------------------------------------------------------------------------- ---------------------------------------------------

There is a rand.c source file in the Redis source code, which is obviously a file related to (pseudo) random numbers. Look at the file code only a few 50 lines, but the principle of the algorithm involved is not simple, although it is somewhat obscure to read, but for in-depth understanding of the 48-bit space pseudo-random number algorithm is a rare model. The author writes in a note to the file that this pseudo-random number generation function is derived from drand48 () in the Pysam source. About Pysam is the project, not the focus, in fact, many UNIX systems exist drand48 this function (svr4,posix.1-2001), we can in the terminal man drand48.

You can see that a series of functions related to random numbers are defined in the header file Stdlib.h, with the exception of Drand48 (), erand48 (), lrand48 (), nrand48 (), and so on. The so-called 48 refers to the bits number of random number seeds. Their main operations are using the same random number generation algorithm, but subsequent processing is different, such as drand48 (), erand48 () returns a floating-point between [0.0, 1.0), while lrand48 () and nrand48 () return [0, 2^31] The integer between. There is an important function--redislrand48 () in the Rand.c file. Although the author is inspired by the drand48 () source code to write this function, but in fact redisLrand48 () and lrand48 () more like, have almost the same programming interface, and set the same random number seed, generated the same pseudo-random number sequence.

Continuing to read the comments at the beginning of the RAND.C, you can see that the author is writing this file to replace the Math.random () function implemented in Lua. The author argues that the libc rand () function is used in the default implementation of the LUA random number generation algorithm, and that libc rand () does not guarantee that the same sequence of random numbers can be generated for different platforms with the same seed set. So the author rewrites this file to ensure cross-platform.

----------------------------------------------------------------------------------------------------------- ---------------------------------------------------

Algorithm principle linear congruence equation the pseudo-random number generation algorithm uses a linear congruence equation, as its name implies: "Linear", refers to the independent variable (x in the equation) and the dependent variable (y in the equation) is a linear relationship; so-called "Congruence", indicating that they are divided by the same number (for example, M) and the remainder is the same. Can be expressed as: Y = f (x)%m or y = f (x) mod M "F (x) is a linear function, that is, X is 1" in Linux, through the LRAND48 manual page, you can see the specific contents of the linear congruence equation used by lrand48 () (Redis uses the same An equation):

Xn + 1 = (aXn + c) mod m, where n> = 0
//by default
a = 0x5DEECE66D
c = 0xB
m = 2 ^ 48
X0, X1, X2 ... Xn is a pseudo-random number sequence. Each of them must be generated by the operation of the last number.
Seeds Random number generation algorithms all have the concept of seeds, whether they are ordinary rand () or lrand48 (). We mentioned the term pseudo-random number many times before. The reason why we say “pseudo” is because after the seed is determined, the sequence of random numbers generated after that is also determined. In other words, as long as it is the same seed, the generated random number sequence is always fixed (under the same platform), not completely random. Through the linear congruence equation in the previous section, we observe that the generation of random number sequences is affected by X0, a, and c (m is always 2 ^ 48, indicating operations in 48-bit space, so these function suffixes will have 48). Therefore, X0, a, c are essentially seeds. These three parameters are provided with the same default values in the implementation of POSIX's lrand48 () and Redis. But the difference is that the POSIX standard provides the lcong48 () function to modify the default values of these three parameters, and Redis can only modify the value of X0 (representing the first number of the random number sequence), so the real There is only one-X0. Three static global variables are defined in the source code implementation variable rand.c: x, a, c correspond to the parameters in the linear congruence equation
static uint32_t x [3] = {X0, X1, X2}, a [3] = {A0, A1, A2}, c = C;
Several macro definitions are involved here:
#define X0 0x330E
#define X1 0xABCD
#define X2 0x1234
#define A0 0xE66D
#define A1 0xDEEC
#define A2 0x5
#define C 0xB
        It can be seen that the three elements in the array a hold the three parts of the constant 0x5DEECE66D in the linear congruence equation. Note that each element in a is of type uint32_t, which means 4 bytes. In fact, we split 0x5DEECE66D into three parts (0x5, 0xDEEC, 0xE66D), and each part occupies at most 2 bytes. Another thing to note is that the semantics of X0 in the macro definition are different from the X0 in the previous section (the previous section). X0 here represents the content of the first element of the array X when it is initialized by default. The received X0 represents the first random number, which is equivalent to X2 * 2 ^ 32 + X1 * 2 ^ 16 + X0 in this section.
External functions rand.c only provides two functions externally:
redisLrand48 (): returns a random number
redisSrand48 (): Set random number seed
There is also a static function next (), which serves redisLrand48 () and is not visible to the outside world. redisSrand48 () function
void redisSrand48 (int32_t seedval) {
    SEED (X0, LOW (seedval), HIGH (seedval));
}
        To set the seed of random numbers, the macro function SEED is used. Look at the definition:
#define SET3 (x, x0, x1, x2) ((x) [0] = (x0), (x) [1] = (x1), (x) (2) = (x2))
#define SEED (x0, x1, x2) (SET3 (x, x0, x1, x2), SET3 (a, A0, A1, A2), c = C)
        It can be seen that the macro function SEED has performed three assignment operations, namely, to the array x, the array a and the variable c. But a and c are assigned default values, so the actual operation of SEED (X0, LOW (seedval), HIGH (seedval)) is to use X0 (0x330E), the lower 16 bits of seedval, and the upper 16 bits of seedval to assign values Give x [0], x [1], x [2]. LOW and HIGH are also two macro functions
#define N 16
#define MASK ((1 << (N-1)) + (1 << (N-1))-1) // MASK = 65535 (16 binary bits 1)
#define LOW (x) ((unsigned) (x) & MASK) // LOW (x) gets the lower 16 bits of 3x (2 bytes)
#define HIGH (x) LOW ((x) >> N) // HIGH (x) gets the upper 16 bits of x (2 bytes)
Note that the upper 16 bits and lower 16 bits refer to the high and low of 32-bit integer data. The final result to be retained each time is 48 bits (stored in the array x []), which is divided into high 16 bits (x [2]), medium 16 bits (x [1]) and low 16 bits (x [ 0]). Don't be confused.

redisLrand48 () function
int32_t redisLrand48 () {
    next ();
    return (((int32_t) x [2] << (N-1)) + (x [1] >> 1));
}
        The next function will be included later, first look at the return value, and bring N = 16 into it, that is:
((int32_t) x [2] << 15) + (x [1] >> 1)
It is x [2] * 2 ^ 15 + x [1] / 2. There is no reason to say here, random numbers, this is just a scheme (the scheme adopted by lrand48 ()). next () function Formula derivation and simplification next () is a static function serving redisLrand48 (), which means that it is not visible to the outside world, but the next () function is the essence of the pseudo-random number generation algorithm. Recall the linear congruence formula again:
Xn + 1 = (aXn + c) mod m, where n> = 0, m = 2 ^ 48
        Note: n and n + 1 are X subscripts. From a mathematical point of view, it is very simple, first multiply and then add, and finally take the modulus. But the computer does not do that, because it will overflow, so we use the array to store x and a. So what you have to do is use arrays to simulate the multiplication of two numbers. First express X and a in the equation:


Convert the linear congruence equation:
Xn + 1 = (aXn + c) mod m
Xn + 1 = ((aXn) mod m + c mod m) mod m
Xn + 1 = ((aXn) mod m + c) mod m
        First calculate a * Xn:
        This polynomial consists of 9 parts. Calculate (a * Xn) mod m again. Because m = 2 ^ 48, in the above polynomial, the term of the power of 2 greater than or equal to 48 must be a multiple of 2 ^ 48, and the modulus is 0, so these terms can be removed. So the result of (a * Xn) mod m is:
        Only 6 items left. Then merge similar items and add c, it is easy to write the result of a * Xn + c:
        The content in each parenthesis of the above formula is the new content saved in x [2], x [1], x [0] after a random number is generated (that is, used to represent Xn + 1). Of course, when programming to calculate, we must also consider the issue of carry. First look at a macro definition: the macro used
#define MUL (x, y, z) {int32_t l = (long) (x) * (long) (y); (z) [0] = LOW (l); (z) [1] = HIGH (l );}
#define CARRY (x, y) ((int32_t) (x) + (long) (y)> MASK)
#define ADDEQU (x, y, z) (z = CARRY (x, (y)), x = LOW (x + (y)))
        Obviously the macro function MUL (x, y, z) represents a multiplication operation, x and y are multiplied, and the result is stored in z. z [0] stores the lower 16 bits, and z [1] stores the upper 16 bits. CARRY (x, y) is to determine whether the addition of two numbers will carry ", here" carry "refers to the size of more than 2 bytes (16 bits). Because we are either array a or array x, each of them The element holds 2 bytes of valid data. That is to say, if the number in x [0] exceeds 2 bytes, the extra part must be "carried" into x [1]. Each The elements are actually 4 bytes in size. The reason why all 4 bytes are not used for storage is because when all 4 bytes are used to store data, the sum may overflow and the compiler will pop up a warning , The result will also become negative.
        The operation performed by ADDEQU (x, y, z) is: add x and y, and store the result in x. Whether to save "carry" in z. next () source code
static void next (void) {
    uint32_t p [2], q [2], r [2], carry0, carry1;

    MUL (a [0], x [0], p);
    ADDEQU (p [0], c, carry0);
    ADDEQU (p [1], carry0, carry1);
    MUL (a [0], x [1], q);
    ADDEQU (p [1], q [0], carry0);
    MUL (a [1], x [0], r);
    x [2] = LOW (carry0 + carry1 + CARRY (p [1], r [0]) + q [1] + r [1] +
            a [0] * x [2] + a [1] * x [1] + a [2] * x [0]);
    x [1] = LOW (p [1] + r [0]);
    x [0] = LOW (p [0]);
}
        Let's analyze it line by line. Do not read the statement, read from line 4.
    MUL (a [0], x [0], p);
    ADDEQU (p [0], c, carry0);
    ADDEQU (p [1], carry0, carry1);
MUL operation: a [0] and x [0] are multiplied, and the result is saved in the array p.
The first ADDEQU: Then p [0] is added to c again, carry0 indicates whether there is a "carry". Here is the operation.
After the second ADDEQU is to add the lower 16 bits to c or not, add to p [1], and determine whether there is a new carry after this addition (refer to the 2 ^ 16 degree polynomial to the 2 ^ 32 degree polynomial Carry) and save to carry1.
At this point, the operations related to the first-degree polynomial are counted, and the next operation is the 2 ^ 16-degree polynomial operation.


    MUL (a [0], x [1], q);
    ADDEQU (p [1], q [0], carry0);
    MUL (a [1], x [0], r);
The first MUL operation: a [0] and x [1] are multiplied, and the result is saved in the array q.
ADDEQU operation: add p [1] (the value carried to the "middle 16 bits" after the "lower 16 bits" operation) and q [0], and save the result to p [1], carry0 identifies the "middle 16 bits" "(2 ^ 16 degree polynomial) whether there is a new carry for the higher 16 bits (2 ^ 32 degree polynomial).
The second MUL operation: a [1] and x [0] are multiplied, and the result is saved in the array r.
So far, the coefficient correlation calculation is completed.
    x [2] = LOW (carry0 + carry1 + CARRY (p [1], r [0]) + q [1] + r [1] +
            a [0] * x [2] + a [1] * x [1] + a [2] * x [0]);
    x [1] = LOW (p [1] + r [0]);
    x [0] = LOW (p [0]);
        Read it backwards, first look at x [0], which should store the result of a first-degree polynomial. x [1] stores the coefficients of the 2 ^ 16 degree polynomial. p [1] stores the result of a [0] x [1] and the lower 16-bit carry. r [0] stores the result of a [1] x [0] (not including the carry).
        x [2] is the most complex. Recall the coefficients of polynomials of degree 2 ^ 32: Let ’s take a look at that line of code, and read the parts that are directly multiplied in pairs, then you can remove them and look at the other parts: carry0 + carry1 + CARRY (p [ 1],
 r [0])
 + q [1]
 + r [1] The remaining polynomial is the carry of the 2 ^ 16 degree polynomial during the operation.
        I drew a picture to show you:
        The upper two rows of this figure are the table headers, and the colored parts below are the stored contents. For example, the third column indicates that x [0] stores the sum of p [0] and c. The transition from the first picture to the second picture is due to the implementation of ADDEQU (p [1],
  q [0],
  carry0); This statement adds q [0] to p [1]. Then look at the code of x [2], it will be much clearer, the point is to consider the carry.

See the pseudo-random number generation algorithm in the Redis source code