本文讨论的情况：已有一个性能符合要求的伪随机数生成器（例如Xorshift128+、MersenneTwister的各种实现），使用该随机数生成器生成随机序列，并映射到更小的取值范围。关于伪随机数生成器（伪随机数算法）本身，以后如果有机会，另开文章来讨论。

0. 随机数之“伪”

“偶然性的一种形式，具有某一概率的事件集合中的各个事件所表现出来的不确定性。”

当我们有一个成熟的伪随机数算法后，其实问题才刚刚开始。本文只讨论最简单的线性分布随机数。

为了方便描述，先做如下定义：

• 对于一个只包含整数的闭区间，定义其包含的整数个数为区间的宽度，例如[0, 25]的宽度为26，[0, 2ⁿ-1]的宽度为2ⁿ，[a, b]的宽度为(b - a + 1)；
• Random()为伪随机数发生器，每次调用Random()都会生成一个指定区间（例如[0, 2ⁿ-1]）内的数，该区间则称为随机值域；
• f(a, b)为一个生成指定区间[a, b]内随机数的函数，则该[a, b]区间则称为目标值域；
• 随机值域中任一个数，都能通过指定动作对应到目标值域中一个数，该动作称为映射；
• 所有元素的理论生成概率中，最大概率与最小概率的比值，称为概率失衡比；

常见的伪随机数算法，其值域一般为[0, 2ⁿ-1]范围内的整数值，对应的二进制表达即n bit的全0~全1的所有值。这个范围对于大多数场景并不直接适用。例如，我希望生成一个全由小写字母组成的随机序列，理论上讲应该生成[0, 25]或者[1, 26]范围的随机数，然后每一个数与一个字母对应。容易想到，对算法生成的随机数对目标值域宽度求模，即可映射到目标值域中的一个数。即：

f(a, b) = Random() mod (b - a + 1) + a

这个方法，实际上会对结果的随机性产生一些瑕疵。设p(x)为x产生的概率，取n = 32, a = 0, b = 25，即以随机数算法值域为[0, 2³²-1]、生成[0, 25]为例：

∵2³² = 26 *  165191049 + 21，

∴p(0)=p(1)=……=p(20)=(165191049 + 1) / 2³² ≈ 0.0384615385,

p(21)=p(22)=……=p(25)=165191049 / 2³² ≈ 0.0384615383

p(0)虽然略大于p(21)，但是这个概率差值很小，可以忽略不计。

考虑一种极端情况，取n = 32, a = 0, b = 2³²-2时：

∵2³² = 1*  (2³²-2)+ 1，

∴p(0)=(1 + 1) / 2³², p(1)=p(2)=……=p(2³²-2)=1 / 2³²,

∴p(0) / p(1) = 2

第一个数的概率是其他数的两倍，这个概率偏差对于伪随机数算法来说，完全无法接受，甚至应该归为错误。

再考虑一种极端情况，取n = 32, a = 0, b = 2³¹时：

∵2³² = 1*  (2³¹+1)+ (2³¹-1)，

∴p(0)=p(1)=……=p(2³¹-1)=(1 + 1) / 2³², p(2³¹)=1 / 2³²

∴p(0) / p(2³¹) = 2

这下更过分，前面的一大堆数的概率都是最后一个的2倍，这个概率偏差对于伪随机数算法来说，完全无法接受，应该归为错误。

更一般地，对于任意取值的n, a, b(a <= b, n > (b - a + 1))，除非n正好能被(b - a + 1)整除，其他情况总有一些数出现的频次比另外一些多1。如果用η表示概率失衡比，⌊a⌋表示实数a的整数部分，则：

η = (⌊n / (b - a + 1)⌋ + 1) / ⌊n / (b - a + 1)⌋ = 1 + 1 / ⌊n / (b - a + 1)⌋

即⌊n / (b - a + 1)⌋越小，偏差越大；也就是说目标值域与随机值域两者的宽度越接近，产生随机序列间的概率失衡越严重。

2. 平衡概率的悖论

从上述推理中可知，对随机值域中的数进行求模，而映射到目标值域的方法，当且仅当随机值域宽度正好可被目标值域宽度整除时，所有元素的出现概率完全一致；否则不可避免将产生概率上的失衡。随机值域与目标值域的关系，可以表示为如下的二维表：

说明如下：

• 表共有m个元素，对应随机值域的宽度；
• 表共有n列，对应目标值域的宽度；
• 表必定包含k个完整行，这k行统称为“体”；
• 当随机值域宽度不能被目标值域宽度整除时，表还包含k+1行，该行为不完整行，称为“尾”；

结合前述的数学计算，可以对应产生以下几条图形特征推论：

• 表是一个完整矩形⟺随机值域宽度可被目标值域宽度整除；
• 表是一个完整矩形⟺所有元素的概率均衡；
• 如果表包含尾，则概率失衡比 = 表行数 / 体行数；

再次考察求模法映射目标值域的过程，可以发现，造成随机概率失衡的原因，正是尾的存在。于是有了以下尝试。

2.1 断尾法

既然惹祸的是尾，那我直接把尾去了，所有元素的概率不就均衡了？程序上倒也不难实现。先计算出体的最大值，只要生成的随机数大于该值则直接舍去，然后直接执行下一次随机数生成。仍然以[0, 2³²-1]为随机值域，则代码实现大概是这个样子：

uint32_t GetRandomUInt32(uint32_t a, uint32_t b)
{
    if(a > b)
    {
        return GetRandomUInt32(b, a);
    }
    if(a == b)
    {
        return a;
    }
    if(0 == a && UINT32_MAX == b)
    {
        return Random();
    }

    uint32_t iRangeSize = b - a + 1;
    uint32_t iTailSize = (UINT_MAX - iRangeSize + 1) % iRangeSize;
    uint32_t iMaxBodyValue = UINT_MAX - iTailSize;

    uint32_t iRet = Random();
    while(iRet > iMaxBodyValue)
    {
        iRet = Random();
    }
    
    return iRet % iRangeSize + a;
}

如此生成，随机性是可靠的，但是有iTailSize / (UINT_MAX + 1)的概率需要重新生成。极端情况下，重新生成的概率接近1/2，相当于需要多耗费将近一倍的资源来生成符合预期的目标值域随机数。

2.2 缩放法

然后又有朋友说，直接将随机值域缩放到目标值域，这样就不会有这个尾了。听上去好像像是那么回事，而且代码实现起来更简单：

uint32_t GetRandomUInt32(uint32_t a, uint32_t b)
{
    if(a > b)
    {
        return GetRandomUInt32(b, a);
    }
    if(a == b)
    {
        return a;
    }
    if(0 == a && UINT32_MAX == b)
    {
        return Random();
    }

    return (float)Random() / UINT_MAX * (b - a) + a;
}

问题在于，随机值域不是连续的。如此缩放后，虽然仍可在数轴上找到对应的点，但是这些点很大一部分将不是整数。此时不论四舍五入或者直接舍弃小数部分，都会导致有一些目标值域中的数被多个随机值域的数映射到，元素的概率仍然是失衡的。

2.3 尾递增法

既然断尾法可能造成资源浪费，那就想办法把尾用回来。比如说：如果随机数落在尾内，使用已知序列来替代。代码实现如下：

uint32_t g_iTailNum = 0;
uint32_t GetRandomUInt32(uint32_t a, uint32_t b)
{
    if(a > b)
    {
        return GetRandomUInt32(b, a);
    }
    if(a == b)
    {
        return a;
    }
    if(0 == a && UINT32_MAX == b)
    {
        return Random();
    }

    uint32_t iRangeSize = b - a + 1;
    uint32_t iTailSize = (UINT_MAX - iRangeSize + 1) % iRangeSize;
    uint32_t iMaxBodyValue = UINT_MAX - iTailSize;

    uint32_t iRet = Random();
    if(iRet > iMaxBodyValue)
    {
        if(UINT_MAX == g_iTailNum)
        {
            iRet = g_iTailNum;
            g_iTailNum = 0;
        }
        else
        {
            iRet = g_iTailNum++;
        }
    }
    
    return iRet % iRangeSize + a;
}

如此生成，只要是随机数生成在尾部，则按顺序依次输出，元素的出现概率被强制保证均衡，但是随机性仍不可靠。极端情况下，可能会有0,1,2,3,4,……这样极不随机的输出。

3. 走马灯的启示

几种对求模法的改进尝试，效果都不尽如人意。好在绝大多数实际使用场景，目标值域的宽度远小于随机值域，这个问题也就搁置了许久，没有深究。转机竟然来自近期偶然看到的一个走马灯非遗节目。节目本身中规中矩，不过走马灯头追尾的循环，给了我一个启示。

如前文所述，随机值域可以表示为一张二维表，表中元素与值域中的整数一一对应。包含尾的表，最末一行一定是不完整的。以生成全小写字母随机序列为例，随机值域的每一个元素与字母的映射关系如下图所示：

其中，绿色部分为随机值域覆盖部分。根据前述图形特征推论，如果绿色部分是一个完整矩形，这就意味着随机值域可以被目标值域整除，概率均衡。对于一张固定的表，只要它存在尾，不论做任何调整，都无法完全抹去最后一行的影响，这个问题似乎无解。如果让这张表动起来，像走马灯一样，头追尾，结果又如何呢？

每一个同色区域都表示某一次调用Random()的映射关系，而此映射关系不是固定的，每一映射关系都是在前一映射关系的基础上偏移指定个元素实现的（头追尾）。通过计算可知，当调用次数为26的整数倍时，所有元素的理论概率达到一次均衡。对于一般情况，只要调用次数足够多，这张由多个动态表组成的大映射表将恢复为标准矩形，那么元素的总体概率将趋向均衡状态。代码实现如下：

uint32_t g_iRangeSize = UINT32_MAX;
uint32_t g_iTailSize = 0;
uint32_t g_iMaxBodyValue = UINT32_MAX;
uint32_t g_iOffset = 0;
uint32_t GetRandomUInt32(uint32_t a, uint32_t b)
{
    if(a > b)
    {
        return GetRandomUInt32(b, a);
    }
    if(a == b)
    {
        return a;
    }
    if(0 == a && UINT32_MAX == b)
    {
        return GetRandomUInt32();
    }

    if(g_iRangeSize != b - a + 1)
    {
        g_iRangeSize = b - a + 1;
        g_iTailSize = (UINT_MAX - g_iRangeSize + 1) % g_iRangeSize;
        g_iMaxBodyValue = UINT_MAX - g_iTailSize;
    }

    uint32_t iRet = a + (Random() - g_iRangeSize + g_iOffset) % g_iRangeSize;

    g_iOffset += g_iTailSize;
    if(g_iOffset >= g_iRangeSize)
    {
        g_iOffset %= g_iRangeSize;
    }

    return iRet;
}

显然，走马灯法对于概率失衡的修正需要通过多次调用逐步实现。极端情况，需要达到极大的调用次数才能实现一次理论概率均衡。本文只是对随机理论的一种讨论，实际使用时，需要综合考虑性能和精度。在绝大多数情况下，断尾法就足够满足使用了。

4. 一点感悟