生成一个不偏的随机整数的最佳算法是什么？该整数需要在一个给定范围内。

Question

在这个 StackOverflow 的问题中： 生成一个指定范围内的随机整数 被接受的答案提供了以下公式，用于生成介于给定的 min 和 max 之间的随机整数，其中min和max被包括在范围内：

output = min + (rand() % (int)(max - min + 1))

但它也说：

这仍然对较小的数字稍微有偏差……也可以扩展它以消除偏差。

但它并没有解释为什么会对较小的数字有偏差，也没有说明如何消除这种偏差。所以问题是：在不依赖于任何花哨的东西，只使用rand()函数生成一个（有符号）范围内的随机整数时，这是最优方法吗？如果是最优的，如何消除偏差呢？

编辑:

我刚刚测试了@Joey建议的while循环算法与浮点外推算法：

static const double s_invRandMax = 1.0/((double)RAND_MAX + 1.0);
return min + (int)(((double)(max + 1 - min))*rand()*s_invRandMax);

为了观察在一定数量的“桶”中，有多少个均匀的“球”正在“落下”并分配，我进行了一个用于浮点数推断和一个用于while循环算法的测试。但结果显示，由于“球”（和“桶”）的数量不同，所以我无法轻易地选择一个赢家。可以在这个 Ideone 页面找到工作代码。例如，使用10个桶和100个球时，浮点数推断的最大偏差比while循环算法小（分别为0.04和0.05），但是当有1000个球时，while循环算法的最大偏差较小（分别为0.024和0.011），而当有10000个球时，浮点数推断再次表现更好（分别为0.0034和0.0053），等等，没有太多的一致性。考虑到可能没有任何一种算法能够比另一种算法更一致地产生更好的均匀分布，这使我更倾向于使用浮点数推断，因为它似乎比while循环算法执行得更快。那么选择浮点数推断算法是否合适，或者我的测试/结论并不完全正确？

- Desmond Hume

---

由于rand()函数生成的数值范围由RAND_MAX定义，而且RAND_MAX通常不能被除数整除，所以会对较小的数字有偏向性。因此，在0到RAND_MAX % divisor - 1之间的所有数字被选中的概率更高。请注意，这里不是解释，只是翻译。 - nhahtdh

假设你有一个随机数生成器，它以相等的概率给出0、1或2。如果你对它应用模运算2，得到0或1，则可以看出0被选中的概率是1的两倍。我想这就是你引用的那个声明所暗示的意思。 - ereOn

http://eternallyconfuzzled.com/arts/jsw_art_rand.aspx - Fred Foo

输入是一个整数，输出也是一个整数 - 在转换过程中转换为浮点数并不会改变某些数字比其他数字更有可能出现的事实，尽管它可能会改变那些数字是哪些。 - Mark Ransom

可能是生成指定范围内的随机整数的重复问题。 - Adrian McCarthy

Answer 1

问题在于你正在进行模运算。如果 RAND_MAX 能够被你的模数整除，那么这就不是一个问题，但通常情况下并非如此。举个非常牵强的例子，假设 RAND_MAX 为11，而你的模数为3。你将得到以下可能的随机数和相应的余数：

0 1 2 3 4 5 6 7 8 9 10
0 1 2 0 1 2 0 1 2 0 1

正如您所看到的，0和1比2稍微更有可能。

解决这个问题的一个选项是拒绝抽样：通过禁止以上的数字9和10，可以使得结果分布再次变为均匀。棘手的部分是如何高效地实现。一个非常好的例子（让我花了两天时间才理解它为什么有效）可以在Java的java.util.Random.nextInt(int)方法中找到。

Java算法的一点棘手之处在于，他们避免使用乘法和除法等缓慢的操作进行检查。如果您不太在意，也可以用朴素的方法来做。

int n = (int)(max - min + 1);
int remainder = RAND_MAX % n;
int x, output;
do {
  x = rand();
  output = x % n;
} while (x >= RAND_MAX - remainder);
return min + output;

编辑：已经修正了上面代码中的一个偏差错误，现在它可以正常工作了。我还创建了一个小样例程序（使用C#；使用一个均匀分布的伪随机数生成器生成0到15之间的数字，并通过不同的方法构建出一个生成0到6之间数字的伪随机数生成器）：

using System;

class Rand {
    static Random r = new Random();

    static int Rand16() {
        return r.Next(16);
    }

    static int Rand7Naive() {
        return Rand16() % 7;
    }

    static int Rand7Float() {
        return (int)(Rand16() / 16.0 * 7);
    }

    // corrected
    static int Rand7RejectionNaive() {
        int n = 7, remainder = 16 % n, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x >= 16 - remainder);
        return output;
    }

    // adapted to fit the constraints of this example
    static int Rand7RejectionJava() {
        int n = 7, x, output;
        do {
            x = Rand16();
            output = x % n;
        } while (x - output + 6 > 15);
        return output;
    }

    static void Test(Func<int> rand, string name) {
        var buckets = new int[7];
        for (int i = 0; i < 10000000; i++) buckets[rand()]++;
        Console.WriteLine(name);
        for (int i = 0; i < 7; i++) Console.WriteLine("{0}\t{1}", i, buckets[i]);
    }

    static void Main() {
        Test(Rand7Naive, "Rand7Naive");
        Test(Rand7Float, "Rand7Float");
        Test(Rand7RejectionNaive, "Rand7RejectionNaive");
    }
}

以下是结果（粘贴到Excel中并添加条件着色以使差异更加明显）：

现在我已经修复了拒绝采样中的错误，它现在按照预期工作（之前会导致偏差为0）。正如你所看到的，浮点数方法并不完美，它只是以不同的方式分布有偏差的数字。

- Joey

---

如果使用(double)rand() / RAND_MAX * n，你会遇到同样的问题，只是更倾向于分布在整个区间上而不是靠近较小的数字。这种方法并不能完全消除偏差。你仍然面临将10个输入数字均匀地分配到3个输出数字中的问题，这是不可能实现的。 - Joey

不要像模数那样将更有偏见的数字放在前面，你可以以0和2，或1和2为例来进行偏置。这有点不可预测，因为它取决于浮点数的外观以及它们何时切换到另一个数字。不过你仍然有偏差。是的，一共有11个数字。对我无法编辑的注释中的小错字请原谅。重点仍然没有改变。 - Joey

@Јοеу 很棒的回答，谢谢。我已经试图理解 nextInt(n) 两天了，直到我偶然发现这个帖子和你的回答。关于你的基准测试的一个评论; 拒绝算法的朴素实现和Java实现之间确实没有太大区别，特别是在除法/乘法操作方面没有更多/更少。或者我漏掉了什么？ - posdef

@posdef：没有。只是Java的变体需要相当多的思考才能理解它的作用，而天真的版本更加直观（我猜应该会稍微慢一些）。 - Joey

@Јοеу 但是由于它们之间实际上没有显着的区别，我不确定为什么您会期望它变慢。 - posdef

Answer 2

问题出现在随机数生成器的输出数量(RAND_MAX+1)不能被期望范围(max-min+1)整除时。由于随机数到输出之间存在一致的映射，某些输出将被映射到比其他输出更多的随机数上。无论如何进行映射 - 可以使用取模、除法、转换为浮点数，或者任何你能想到的方法，基本问题仍然存在。

问题的规模非常小，不苛求的应用程序通常可以忽略它。范围越小，RAND_MAX越大，影响就越不明显。

我拿了你的示例程序并进行了一些调整。首先，我创建了一个特殊版本的rand，只有0-255的范围，以更好地演示效果。我对rangeRandomAlg2进行了一些调整。最后，我将“球”的数量改为1000000以提高一致性。你可以在这里看到结果：http://ideone.com/4P4HY 请注意，浮点版本会产生两个紧密分组的概率，接近0.101或0.097，中间没有其他选项。这就是偏差在起作用。

我认为称其为“Java算法”有点误导-我确信它比Java要古老得多。

int rangeRandomAlg2 (int min, int max)
{
    int n = max - min + 1;
    int remainder = RAND_MAX % n;
    int x;
    do
    {
        x = rand();
    } while (x >= RAND_MAX - remainder);
    return min + x % n;
}

- Mark Ransom

---

嗯，当我阅读Java源代码时，我发现了他们的实现（尽管有点不太明显，但我仍然觉得相当不错）。当然，这并不是任何答案中所写的，因为你从我的答案中复制的那个天真的检查要容易理解得多。 - Joey

@AmbrozBizjak，最坏情况是当n等于RAND_MAX/2+1时，平均会丢弃一半的随机数。当然，最坏情况可能会更糟，但如果运行足够多次，它应该回到平均水平。如果一个随机数生成器没有在其范围内产生所有可能的值，它将无法通过大多数随机性定义的测试。 - Mark Ransom

最坏情况下，n次循环的平均次数恰好为2，因为它是概率1 + 1/2 + 1/4 + 1/8 + 1/16 + 1/32 + ...的无限和。然而，这只是循环次数的平均值。实际的循环次数是不确定的。你可能会真的循环50次，但这是极其罕见的。 - Todd Lehman

@DesmondHume 当我第一次回答这个问题时，我还没有看过《迷失》。你的用户名现在更有趣了。 - Mark Ransom

@MarkRansom，rand() 可能的输出数量等于 RAND_MAX + 1，即 rand() 生成的伪随机数在范围 0 到 RAND_MAX（包括两端）之间。因此，当 min == 0 且 max == RAND_MAX 时，在一些系统中会发生有符号整数溢出（当 RAND_MAX == INT_MAX 时），导致未定义行为，尽管参数 min == 0 和 max == RAND_MAX 是完全有效的。在 RAND_MAX < INT_MAX 的系统中，该函数仍然会进入无限循环。 - Kushagr Jaiswal

Answer 3

很容易看出为什么这个算法会产生偏差的样本。假设你的rand()函数返回集合{0,1,2,3,4}中的均匀整数。如果我想用它来生成一个随机位0或1，我会使用rand() % 2。集合{0,2,4}给我0，而集合{1,3}给我1 - 所以很明显我用60%的概率采样0和40%的概率采样1，根本不是均匀的！

要解决这个问题，你必须确保你所需的范围可以被随机数生成器的范围整除，否则就要在随机数生成器返回大于目标范围的最大可能倍数的数字时丢弃结果。

在上面的例子中，目标范围是2，适合随机生成范围的最大倍数是4，因此我们丢弃任何不在集合{0,1,2,3}中的样本，并重新生成。

- Kerrek SB

Answer 4

到目前为止，最简单的解决方案是使用 std::uniform_int_distribution<int>(min, max)。

- MSalters

Answer 5

你提到了一个关于随机整数算法的两个问题：它是否最优，以及它是否无偏？ 最优 有许多方法可以定义“最优”算法。在这里，我们从平均使用的随机位数来看“最优”算法。在这个意义上，rand不适合用于随机生成数字，因为除了rand()的其他问题>, 它不一定会产生随机位（因为RAND_MAX没有完全指定）。相反，我们将假设我们有一个“真正”的随机生成器，它可以产生无偏和独立的随机位。

1976年，D. E. Knuth和A.C. Yao证明了任何只使用随机位生成给定概率的随机整数的算法都可以表示为二叉树，其中随机位指示遍历树的方向，每个叶子节点（终点）对应一个结果。（Knuth和Yao，“The complexity of nonuniform random number generation”，收录于《Algorithms and Complexity》，1976年。）他们还给出了一个给定算法平均需要的比特数的上下界。在这种情况下，生成[0，n)范围内均匀分布的整数的最优算法将平均需要至少log2(n)个比特，最多log2(n)+2个比特。

在这个意义下，有许多最优算法的例子。请参见我的以下回答：

• 如何从一串随机位中生成[0，n]范围内的随机整数而不浪费位？

无偏

然而，任何同时是最优整数生成器和无偏的生成器，在最坏情况下通常会无限运行，正如Knuth和Yao所示。回到二叉树的问题上，每个n个结果都标记了二叉树中的叶子节点，以便[0，n)中的每个整数都可以以1/n的概率出现。但是，如果1/n有一个非终止的二进制扩展（如果n不是2的幂，则这将是情况），则这棵二叉树必须要么具有“无限”的深度，要么包括树末端的“拒绝”叶子节点。在任一种情况下，该算法都无法在恒定时间内运行，并且在最坏情况下将永远运行。（另一方面，当n是2的幂时，最优二叉树将具有有限深度且不包含拒绝节点。）

对于一般的n，没有办法在不引入偏差的情况下解决最坏时间复杂度问题。例如，模数约简（包括你问题中的min + (rand() % (int)(max - min + 1))）等价于一个二叉树，其中拒绝节点被替换为标记结果 - 但由于可能的结果比拒绝节点更多，只有一些结果可以取代拒绝节点，从而引入偏差。如果在一定迭代次数后停止拒绝，将得到相同类型的二叉树 - 以及相同类型的偏差。 （然而，这种偏差可能在应用程序中可以忽略不计。随机整数生成也有安全方面的考虑，这些太复杂了，无法在此回答中讨论。）

- Peter O.

Answer 6

不失一般性，生成[a，b]范围内的随机整数问题可以简化为在[0，s)范围内生成随机整数的问题。从均匀PRNG生成有界范围上的随机整数的现代技术由以下最近的出版物代表：

Daniel Lemire，“Fast Random Integer Generation in an Interval.” ACM Trans. Model. Comput. Simul. 29, 1, Article 3 (January 2019) (ArXiv draft)

Lemire表明他的算法提供无偏结果，并受到非常快速且高质量的PRNG（例如Melissa O'Neill的PCG发生器）越来越受欢迎的推动，展示了如何快速计算结果，几乎总是避免慢速分割操作。

下面展示了他算法的ISO-C实现，可以在randint()中看到。这里我结合了George Marsaglia的旧版KISS64伪随机数生成器进行演示。出于性能原因，所需的64×64→128位无符号乘法通常最好通过特定于机器的内部函数或内联汇编来实现，直接映射到适当的硬件指令。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

/* PRNG state */
typedef struct Prng_T *Prng_T;
/* Returns uniformly distributed integers in [0, 2**64-1] */
uint64_t random64 (Prng_T);
/* Multiplies two 64-bit factors into a 128-bit product */
void umul64wide (uint64_t, uint64_t, uint64_t *, uint64_t *);

/* Generate in bias-free manner a random integer in [0, s) with Lemire's fast
   algorithm that uses integer division only rarely. s must be in [0, 2**64-1].

   Daniel Lemire, "Fast Random Integer Generation in an Interval," ACM Trans.
   Model. Comput. Simul. 29, 1, Article 3 (January 2019)
*/
uint64_t randint (Prng_T prng, uint64_t s) 
{
    uint64_t x, h, l, t;
    x = random64 (prng);
    umul64wide (x, s, &h, &l);
    if (l < s) {
        t = (0 - s) % s;
        while (l < t) {
            x = random64 (prng);
            umul64wide (x, s, &h, &l);
        }
    }
    return h;
}

#define X86_INLINE_ASM (0)

/* Multiply two 64-bit unsigned integers into a 128 bit unsined product. Return
   the least significant 64 bist of the product to the location pointed to by
   lo, and the most signfiicant 64 bits of the product to the location pointed
   to by hi.
*/
void umul64wide (uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo)
{
#if X86_INLINE_ASM
    uint64_t l, h;
    __asm__ (
        "movq  %2, %%rax;\n\t"  // rax = a
        "mulq  %3;\n\t"         // rdx:rax = a * b
        "movq  %%rax, %0;\n\t"  // l = (a * b)<31:0>
        "movq  %%rdx, %1;\n\t"  // h = (a * b)<63:32>
        : "=r"(l), "=r"(h)
        : "r"(a), "r"(b)
        : "%rax", "%rdx");
    *lo = l;
    *hi = h;
#else // X86_INLINE_ASM
    uint64_t a_lo = (uint64_t)(uint32_t)a;
    uint64_t a_hi = a >> 32;
    uint64_t b_lo = (uint64_t)(uint32_t)b;
    uint64_t b_hi = b >> 32;

    uint64_t p0 = a_lo * b_lo;
    uint64_t p1 = a_lo * b_hi;
    uint64_t p2 = a_hi * b_lo;
    uint64_t p3 = a_hi * b_hi;

    uint32_t cy = (uint32_t)(((p0 >> 32) + (uint32_t)p1 + (uint32_t)p2) >> 32);

    *lo = p0 + (p1 << 32) + (p2 << 32);
    *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
#endif // X86_INLINE_ASM
}

/* George Marsaglia's KISS64 generator, posted to comp.lang.c on 28 Feb 2009
   https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
*/
struct Prng_T {
    uint64_t x, c, y, z, t;
};

struct Prng_T kiss64 = {1234567890987654321ULL, 123456123456123456ULL,
                        362436362436362436ULL, 1066149217761810ULL, 0ULL};

/* KISS64 state equations */
#define MWC64 (kiss64->t = (kiss64->x << 58) + kiss64->c,            \
               kiss64->c = (kiss64->x >> 6), kiss64->x += kiss64->t, \
               kiss64->c += (kiss64->x < kiss64->t), kiss64->x)
#define XSH64 (kiss64->y ^= (kiss64->y << 13), kiss64->y ^= (kiss64->y >> 17), \
               kiss64->y ^= (kiss64->y << 43))
#define CNG64 (kiss64->z = 6906969069ULL * kiss64->z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
uint64_t random64 (Prng_T kiss64)
{
    return KISS64;
}

int main (void)
{
    int i;
    Prng_T state = &kiss64;

    for (i = 0; i < 1000; i++) {
        printf ("%llu\n", randint (state, 10));
    }
    return EXIT_SUCCESS;
}

- njuffa

Answer 7

如果你真的想要一个完美的生成器，假设你有一个完美的rand()函数，你需要应用下面解释的方法。

我们将创建一个从0到max-min=b-1的随机数r，然后很容易将其移动到您想要的范围，只需取r+min即可。

我们将创建一个随机数，其中b < RAND_MAX，但该过程可以轻松地采用任何基数的随机数

过程：

1. 在其原始RAND_MAX大小中获取随机数r，不进行任何截断
2. 以基数b显示此数字
3. 对于从0到b-1的m个随机数，取此数字的前m=floor(log_b(RAND_MAX))个数字
4. 将每个数字都移位min（即r+min），以使它们进入所需的范围（min，max）

由于log_b(RAND_MAX)不一定是整数，表示中的最后一位数字被浪费了。

仅使用mod（％）的原始方法是错误的，正是因为

(log_b(RAND_MAX) - floor(log_b(RAND_MAX)))/ceil(log_b(RAND_MAX)) 

你可能会认为这并不多，但如果你坚持要精确，那就是这个过程。

- user6110389