生成n套数位加减乘除_伪随机数生成：梅森旋转（Mersenne Twister/MT）算法

(给算法爱好者加星标，修炼编程内功)

作者：tik_tokc97

https://blog.csdn.net/tick_tock97/article/details/78657851

前言

最近在看吴军博士的《数学之美》一书，把很多之前没注意到，没用到，甚至不知道怎么用的数学知识和实际问题联系了起来，感觉打开了新世界的大门一样。这本书很多知识点还有技术都是点到为止，并没有深入，所谓师傅领进门，修行在个人吧。

在本书第16章《信息指纹及其应用》一文中，介绍到了现在常用的一个伪随机数生成算法，梅森旋转(Mersenne Twister/MT)算法，它的周期长，对于一个k位2进制数，梅森旋转算法可在[0,2^k-1]的范围内生成离散型均匀分布的随机数。

以下内容大多翻译自Wikipedia(中文wiki里的和英文的差的实在太多了……不然我就直接搬运了)

正文

优点

·许可免费，而且对所有它的变体专利免费(除CryptMT外)
·几乎无处不在：它被包含在大多数编程语言和库中
·通过了包括Diehard测试在内的大多数统计随机性测试(除TestU01测试外)
·在应用最广泛的MT19937变体中，周期长达2^19937-1
·在MT19937-32的情况下对1 ≤ k ≤ 623，满足k-分布
·比其他大多数随机数发生算法要快

k-分布

一个周期为P的w位整数的随机序列xi，当满足如下条件时被称为满足v位的k-分布：

假设truncv(x)表示x的前v位形成的数字，并且长度为P的kv位序列：
其中每个可能出现的2^kv组合在一个周期内出现相同的次数(除全0序列出现次数次数比其他序列少1次)

缺点

·需要大量的缓冲器(2.5kib)，但在TinyMT版本中修正
·吞吐量中等，但在SFMT版本中修正
·产生的随机数与seed相关，不能用于蒙特卡洛模拟
·由相同的初始序列产生的随机状态几乎相同，但在2002年的更新中对MT算法的初始化进行了改进，使得对于相同的初始序列也能产生不同的随机状态
·非加密安全的，除CryptMT外

算法详细

本算法基于标准(线性)旋转反馈移位寄存器(twisted generalised feedback shift register/TGFSR)产生随机数

算法中用到的变量如下所示：
    ·w：长度(以bit为单位)
    ·n：递归长度
    ·m：周期参数，用作第三阶段的偏移量
    ·r：低位掩码/低位要提取的位数
    ·a：旋转矩阵的参数
    ·f：初始化梅森旋转链所需参数
    ·b,c：TGFSR的掩码
    ·s,t：TGFSR的位移量
    ·u,d,l：额外梅森旋转所需的掩码和位移量
MT19937-32的参数列表如下：
    ·(w, n, m, r) = (32, 624, 397, 31)
    ·a = 9908B0DF(16)
    ·f = 1812433253
    ·(u, d) = (11, FFFFFFFF16)
    ·(s, b) = (7, 9D2C568016)
    ·(t, c) = (15, EFC6000016)
    ·l = 18
MT19937-64的参数列表如下：
    ·(w, n, m, r) = (64, 312, 156, 31)
    ·a = B5026F5AA96619E9(16)
    ·f = 6364136223846793005
    ·(u, d) = (29, 555555555555555516)
    ·(s, b) = (17, 71D67FFFEDA6000016)
    ·(t, c) = (37, FFF7EEE00000000016)
    ·l = 43

整个算法分为三个阶段(如图所示)：
第一阶段：初始化，获得基础的梅森旋转链；
第二阶段：对于旋转链进行旋转算法；
第三阶段：对于旋转算法所得的结果进行处理；

初始化

首先将传入的seed赋给MT[0]作为初值，然后根据递推式：MT[i] = f × (MT[i-1] ⊕ (MT[i-1] >> (w-2))) + i递推求出梅森旋转链。伪代码如下：

 // 由一个seed初始化随机数产生器 function seed_mt(int seed) {     index := n     MT[0] := seed     for i from 1 to (n - 1) {         MT[i] := lowest w bits of (f * (MT[i-1] xor (MT[i-1] >> (w-2))) + i)     } }

对旋转链执行旋转算法

遍历旋转链，对每个MT[i]，根据递推式：MT[i] = MT[i+m]⊕((upper_mask(MT[i]) || lower_mask(MT[i+1]))A)进行旋转链处理。

其中，“||”代表连接的意思，即组合MT[i]的高 w-r 位和MT[i+1]的低 r 位，设组合后的数字为x，则xA的运算规则为(x0是最低位)：

伪代码为：

lower_mask = (1 << r) - 1upper_mask = ！lower_mask // 旋转算法处理旋转链 function twist() {    for i from 0 to (n-1) {        int x := (MT[i] & upper_mask)+ (MT[(i+1) mod n] & lower_mask)        int xA := x >> 1        if (x mod 2) != 0 {          // 最低位是1            xA := xA xor a        }        MT[i] := MT[(i + m) mod n] xor xA    }    index := 0}

对旋转算法所得结果进行处理

设x是当前序列的下一个值，y是一个临时中间变量，z是算法的返回值。则处理过程如下：
        y := x ⊕ ((x >> u) & d)
        y := y ⊕ ((y << s) & b)
        y := y ⊕ ((y << t) & c)
        z := y ⊕ (y >> l)
伪代码如下：

// 从MT[index]中提取出一个经过处理的值// 每输出n个数字要执行一次旋转算法，以保证随机性 function extract_number() {     if index >= n {         if index > n {           error "发生器尚未初始化"         }         twist()     }     int x := MT[index]     y := x xor ((x >> u) and d)     y := y xor ((y <s)      y := y xor ((y <and c)     z := y xor (y >> l)     index := index + 1     return lowest w bits of (z) }

MT-19937-32实现代码(C语言版)

#include // 定义MT19937-32的常数enum{    // 假定 W = 32 (此项省略)    N = 624,    M = 397,    R = 31,    A = 0x9908B0DF,    F = 1812433253,    U = 11,    // 假定 D = 0xFFFFFFFF (此项省略)    S = 7,    B = 0x9D2C5680,    T = 15,    C = 0xEFC60000,    L = 18,    MASK_LOWER = (1ull <    MASK_UPPER = (1ull <};static  uint32_t  mt[N];static  uint16_t  index;// 根据给定的seed初始化旋转链void  Initialize(const uint32_t  seed){    uint32_t  i;    mt[0] = seed;    for ( i = 1; i     {        mt[i] = (F * (mt[i - 1] ^ (mt[i - 1] >> 30)) + i);    }    index = N;}static void  Twist(){    uint32_t  i, x, xA;    for ( i = 0; i     {        x = (mt[i] & MASK_UPPER) + (mt[(i + 1) % N] & MASK_LOWER);        xA = x >> 1;        if ( x & 0x1 )        {            xA ^= A;        }        mt[i] = mt[(i + M) % N] ^ xA;    }    index = 0;}// 产生一个32位随机数uint32_t  ExtractU32(){         uint32_t  y;    int       i = index;    if ( index >= N )    {        Twist();        i = index;    }    y = mt[i];    index = i + 1;    y ^= (y >> U);    y ^= (y <    y ^= (y <    y ^= (y >> L);        return y;}

- EOF -

推荐阅读   点击标题可跳转

1、随机算法之水塘抽样算法

2、详解各种随机算法

3、随机森林(Random Forest)

觉得本文有帮助？请分享给更多人

推荐关注「算法爱好者」，修炼编程内功

点赞和在看就是最大的支持❤️